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Modelling of green water ingress 
into holds of an open-top containership 
in its preliminary design phase 


Tomasz Cepowski, Assoc. Prof. 
Szczecin Maritime University 


ABSTRACT 


In this paper a method is presented of modelling the green water ingress into holds of open- 
top containership, which can be useful in the preliminary ship design phase. As a result 
of the research a mathematical formula which makes it possible to determine a minimum 
y freeboard height with a view ofas- low- as- possible occurrence rate of green water ingress 


Y 


` into holds at given ship design parameters, was obtained. The research was carried out under 
£ assumption of constant ship hull dimensions. The design formula was elaborated by using 


a method based on a goal-oriented conceptual approach to formulation of design criteria, 
proposed by IMO. On the basis of the concept a deterministic scenario describing operational conditions 
of the ship in question, was assumed, and for the conditions the research was performed. 


Keywords: green water ingress into ship holds, goal-oriented stability standards, modelling, 
open-top containership, preliminary designing, ship hull 


INTRODUCTION 


The traditional design process of ships develops in an 
iterative way in accordance with the so-called Evans or 
Andrews-Evans spiral where design characteristics are 
evolutionarily improved by making successive trial — error 
steps. Knowledge about a designed ship in preliminary design 
phase is slight as several technical ship parameters depend on 
data which are not available in this phase of designing. In order 
to increase degree of accuracy of the preliminary design phase 
the designers must be fitted with appropriate tools for precise 
prediction of selected design characteristics in this design 
phase, based on ship’s main design parameters. 

One ofthe main design problems of open-top containerships 
is green water entering the deck resulting in water ingress 
into the holds. As this type of ships does not comply with 
the requirements of the Load Line Convention, Maritime 
Administration may admit the ship to service at sea provided that 
aditional model tests described in [6] would be performed (the 
test conditions are described in a further part of the paper). 

Amount of water ingress into the holds of open-top 
containership depends a.o. on: wave parameters, ship motion 
parameters, ship hydromechanical and geometrical parameters 
such as: form and dimensions of underwater part of ship’s hull 
(especially ship length), draught, form of ship’s hull stern and 
bow part, freeboard height, deck shape as well as bulwark 
structure [2, 3, 13]. 

Number of geometrical parameters which influence the 
phenomenon of water ingress into the holds is large, hence in 
order to eliminate unimportant ones the investigations can be 
carried out in two ways (according to [11]): 


* to change hull dimensions - keeping hull form unchanged 
(Fig. la) 

* to change hull form - keeping hull dimensions unchanged 
(Fig. 1b). 


a)r] 


w a3 
1 


e a 


Fig. 1. Two ways of investigation of influence of hull form 
and hull dimensions on seakeeping qualities [11]. 
a) unchanged hull form, changeable hull dimensions; 
b) unchanged hull dimensions, changeable hull form 
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Therefore the investigations in question were split in two 

stages: 

ve the investigations of influence of ship hull form parameters 
on green water ingress into the holds at ship hull dimensions 
kept unchanged 

ve the investigations of influence of ship hull dimensions 
on green water ingress into the holds at hull form kept 
unchanged 


In this paper are presented results concerning the first stage 
of the investigations, i.e. the hull form influence on water 
ingress into holds at hull dimensions kept unchanged. 

The investigations were focused only on: 
= influence of form parameters of ship’s hull underwater 

part 
= influence of freeboard height (FB) 
= influence of initial transverse metacentric height (GM) 


out of all parameters influencing the phenomenon of water 
ingress into the holds. 


In compliance with [2, 13] these are the factors most 
influencing the phenomenon of green water entering the 
deck, which - in consequence - are capable of influencing 
water ingress into the holds, to a large extent. And, the initial 
transverse metacentric height constitutes a parameter which 
significantly influences ship rolling, that can in turn influence 
relative ship motion in the deck- side zone, hence also green 
water entering the deck and holds. 

The investigations were performed by applying numerical 
calculations based on the methods not adjusted to determination 
of amount of water ingress into the holds (but only of statistical 
frequency of water entering the deck). From results of the 
research [1, 5] implies the relations between probability of 
water entering the deck and amount of water ingress into the 
holds (Tab. 1). And, the results concerning FPSO ships show 
the linear relation between the water height on the deck and 
freeboard height exceedance (Fig. 2) [2, 3]: 


HO =a: h (1) 
where: 
HO — water height on the deck 
a, — a parameter dependent on location of a given point on 

the deck, and on bow form 
h — maximum exceedance of freeboard [3]: 
h=r—FB (2) 

where: 
r — relative wave motion (Fig. 3) 


FB — freeboard. 


Tab. 1. Water ingress into the holds, where: V „— hourly volumetric rate of 
water ingress, p ,,,— occurrence probability ofh water entering the deck [1, 5] 


Test no. 


V „ [mh] 


Therefore in the presented investigations the linear 
relationship between occurrence rate of green water entering 
the deck and amount of water on the deck, was assumed, and 
the phenomenon of water ingress into the holds was modeled 
by applying the statistical occurrence rate of green water ingress 
onto the deck. 
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Water height on deck (H0) [m] 


0 5 10 15 
Exceedance of freeboard (h) [m] 


Fig. 2. Relation between the water height on the deck, HO, 
and the exceedance of freeboard, h, for a FPSO ship [3] 


: Relative 
i wave 
motion (r) 


Freeboard (fb) | 


Fig. 3. Definition of the relative wave motion (r) and freeboard (fb) [3] 


Hence the investigations presented in this paper were aimed 

at the following: 

> determination of how far initial transverse metacentric 
height influences water ingress into the holds 

> determination of which hull form parameters influence 
significantly water ingress into the holds at the most 
unfavourable values of initial transverse metacentric height 
kept constant 

> determination of a relation between important hull form 
parameters, freeboard height and water ingress into the 
holds 

> elaboration of an analytical function which would make it 
possible to calculate a minimum freeboard height depending 
on hull form parameters and an assumed occurrence rate of 
water ingress into the holds. 


For the investigations the method shown in Fig. 4 was 
chosen. In the proposed method it was assumed that the criterion 
of water ingress into the holds was formulated for the operational 
conditions complying with [9, 11], assuming an appropriate 
deterministic scenario. Such approach is in accordance with 
the conceptual ,,goal-oriented” approach to formulation design 
criteria, proposed by the IMO Maritime Safety Committee [9]. 
As the deterministic scenario (operational conditions) was 
assumed the following conditions which comply with the IMO 
recommendations concerning execution of model tests on open- 
top containerships [6] were selected: 
© the recommended significant wave height H, = 8.5 m 
> wave periods taken every 30 s within the range from 3.5 s 

to 17s 
> the ship motion parameters most unfavourable with a view 

of green water entering the deck: (the ship speed v and wave 

encounter angle f) 
> full load conditions (i.e. maximum ship draught) 
> the investigations were performed for the point on the deck, 

in which the greatest green water ingress onto the deck 

occurs (bow part in vicinity of 1st hold). 


Preliminary numerical calculations were performed to 
determine the ship motion parameters most unfavourable 


with a view of water entering the deck. From the calculations 
resulted that for the assumed ship hull forms (described in 
a further part of this paper) the most unfavourable ship motion 
parameters were the following: the ship speed v = 10 m/s and 
the wave encounter angle B = 60° (in the reference frame: 
0° — following wave, 180° —head wave). As a matter of fact the 
value of the wave encounter angle B is not in compliance with 
the research [1], but may be that for such B values the relative 
motions of underwater part of ship hull are the greatest and 
cause the largest amount of water to enter the deck. 

Ship hull forms obtained by modelling the hull form with 
the use of the curve of frame section areas, were assumed the 
model variants. 


List of model variants 
of ship hull form 


Operational 
conditions 


Design criterion 
of minimum condition 
of water ingress 


Calculation of model values 
of water ingress 
into the holds by using 


accurate methods into the holds 


Analysis of the results 
and determination of crucial 
parameters which affect 
water enetering the holds 


Approximation of model values 
of water ingress into the holds 

with a view of crucial hull form 
parameters and freeboard height 


Elaboration of a relation 
between freeboard 
height and hull 
form parameters 


Fig. 4. Schematic diagram of the investigation method 


MODELLING SHIP’S HULL FORM 
BY MEANS OF THE CURVE OF FRAME 
SECTION AREAS 


In the presented investigations ship’s hull form was 
represented by using the curve of frame section areas F(x). The 
curve was modeled in the same way as presented in [7], i.e. by 
means of the following analytical functions: 


F, =A + Bx + Cx? + Dx? (3) 
F; = E + Fx + Gx? + Hx’ (4) 
where: 
F, — values of the curve of frame section 
areas of stern part of hull (Fig. 5) 
F — values of the curve of frame section 
areas of bow part of hull 


x — distance 
A, B, C, D, E, F, G, H — constants determined under the 
following assumptions (acc. [7]): 


F, =0forx=0 (5) 


dF /dx = AS for x = 0 (6) 
dF ,/dx = 0 forx=L, (7) 
F =F forx=L, (8) 
dF /dx =0 forx=L,,, (9) 
Fe = Faa fOr X = Lam (10) 
FF, =Uforx=L (11) 
dF /dx = FS for x = L (12) 
where: 
L — ship length 
F_ — area of frame section in the midship body zone 


L, — length of stern part of hull 
L,,, — length of stern part and midship body of hull 


Al 
AS. — slope of the curve F(x) at aft perpendicular 
FS -— slope of the curve F(x) at fore perpendicular. 


Lam 


F [m’] 


x [m] 
Fig. 5. The curve of frame section areas F(x) in a parametric form 


The solution of the Eqs. (3) and (4) together with the 
conditions: (5), (6), (7), (8), (9), (10), (11) and (12) made it 
possible to determine the constants as follows: 


A=0 (13) 

B=AS (14) 
D=C2F 7 AS LVL? (15) 
C==AS/2L,)=1,5 DL, (16) 


H=[2F_,.-SF(L-L,,)V(L-L,,,)? 17) 

G= [+ -H@’+2L,,3-3L,,7L)M(L-L,,,? (18) 
F=-2GL u 3H Lyi (19) 
=-FL-GL’?-HL? (20) 


On the basis of the parametrized curve of frame section areas 
the following ship hull form parameters were formulated: 
> length of stern part of hull: L, 
length of stern part and midship body of hull: L y 
slope of the curve F(x) at aft perpendicular: AS 
slope of the curve F(x) at fore perpendicular: FS 
block coefficient of underwater stern part of hull: CBA 


DOG 


e 


KS 


e 


te 


e 


ko 


e 


te 


(21) 


where: 

L — ship length 

V — volumetric displacement of underwater part of hull 
F(x) — parametrized curve of frame section areas 
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“* block coefficient of underwater bow part of hull - CBF: 
I 


[Foodx 


CBE = (22) 
V 


LIST OF MODEL VARIANTS 
OF SHIP’S HULL 


The list of model variants of ship’s hull was elaborated 
under the following limitations: 
+ ship length: L = 144 m 
+ frame section area in the midship body zone: 

max = 173 m? 

+ design draught: d = 8.255 m 
+ length of stern part of hull: L, = 46 + 60 m 
+ length of stern part and midship body of hull: 

L, u7 64 = 95 m 
+ slope of the curve F(x) at aft perpendicular: AS =1+5m 
+ slope of the curve F(x) at fore perpendicular: FS = 1 + 5 m 
+ freeboard heights: FB = 0.75 + 7.75 m. 


Hull block coefficient values were calculated on the basis 
of the above given parameters (in compliance with Eqs. (21) 
and (22)) and amounted to: 
A block coefficient of underwater stern part of hull: 

CBA = 0.53 = 0.61 
A block coefficient of underwater bow part of hull: 

CBF = 0.54 = 0.62. 


On the basis of the above given limitations and Eqs. (3) 
and (4) 16 model variants of underwater part of hull form were 
elaborated. For each of the variants, 5 variants of abovewater 
part of hull form were designed by applying various values of 
freeboard height. As a result, 80 model variants of ship’s hull 
were obtained. 

In the first phase of the research influence of initial 
transverse metacentric height on frequency of water ingress 
into the holds was analyzed. 


INFLUENCE OF INITIAL TRANSVERSE 
METACENTRIC HEIGHT ON 
OCCURRENCE RATE OF WATER 
INGRESS INTO THE HOLDS 


In order to determine influence of initial transverse 
metacentric height on frequency of water ingress into the 
holds calculations were performed by assuming ship motion 
parameters and wave parameters in compliance with the 
assumed deterministic scenario. 

The calcultaions were carried out by means of SEAWAY 
software. This is a computer program based on the theory of 
two-dimensional flow around a body, by which ship motions 
on regular wave and irregular one can be calculated. Its 
validation tests, presented in [8], demonstrate a high accuracy 
of calculations. 

The influence of the initial transverse metacentric height 
GM on occurrence rate of water ingress into the holds was more 
or less similar for all the variants of underwater part of hull 
form and depended first of all on the height of freeboard FB. 
In Fig. 6 the influence of GM on frequency of water ingress into 
the holds is presented for one selected variant of underwater 
part of hull form. 

As results from Fig. 6 the metacentric height GM influences 
ignificantly the water ingress into the hold in the range of 
the freeboard (FB) value equal to or greater than 3.75 m (for 
smaller FB values the influence of GM is negligible). In such 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 
GM [m] 
—FB=7.75m e@FB=3.75m =FB=0.75m 
—FB=5.75m ~ FB=1.75m 


Fig. 6. The influence of the initial transverse metacentric height GM 
on occurrence rate of water ingress into the hold no. 1, FB = var, 
the hull form parameters: Liy = 95m, L, 60 m, AS = 1, FS = 2.2, 
CBA = 0.53, CBF = 0.62, ship speed v = 10 m/s, wave encounter angle 
p = 60°, significant wave height H, = 8.5 m 


cases for GM value greater than about | m significant drop of 
intensity of water ingress into the hold is observed, and the 
greatest intensity occurs when GM = 0.6 m. 

As initial transverse metacentric height is an operational 
parameter dependent a.o. on ship weight and location of 
its centre of gravity, it is hard to consider it to be a design 
parameter. Hence in the further part of the research it was 
decided not to take the parameter into account and to conduct 
calculations of water ingress into the hold of model ship 
variants for a GM value at which the water ingress into the hold 
is the greatest. Such approach allows to decrease number of 
independent variables in regression model and leads to errors 
on the so called „safe side”. 


APPLICATION OF STATISTICAL 
METHODS FOR APPROXIMATION 
OF WATER INGRESS INTO THE HOLDS 


In this part of the research calculations of water ingress 
into the holds were performed for all design ship variants and 
the assumed value of the initial transverse metacentric height 
GM = 0.6 m. The calculations were performed with the use of 
accurate methods (namely the SEAWAY software) and input 
data selected from operational ones resulting from the assumed 
deterministic scenario. 

By making use of statistical methods (namely analysis of 
variance) and those based on artificial neural networks (namely 
sensitivity analysis) it was determined which design parameters 
affect occurrence rate of water ingress into the holds. Next, 
by making use of linear regression method the following 
analytical relation between the crucial design parameters and 
the occurrence rate of water ingress into the holds, nz, was 
determined: 

nz =-  FB-a,:L,+ 


- a, * CBF — o, CBA 


(23) 


where: 

nz -— statistical hourly occurrence rate of water ingress into 
the holds 

FB — freeboard height [m] 

L, — length of ship hull stern part (acc. Fig. 6) [m] 

CBF — block coefficient of underwater ship hull bow part [-] 

CBA — block coefficient of underwater ship hull stern part [-] 


a, = 1111.63 
a, = 25.93 
a, = 5.23 

a, = 229.68 
a, = 729.82. 


The relation (23) is characterized of: 

* high value of the determination coefficient R? = 0.96 
(explained variance) 

¢ high significance of free terms and all independent variables, 
being on the significance level a = 0.05, 

which proves that the proposed relations are properly adjusted 

to the model data. 


The regression of the dependent variable nz indicate 
that all components of Eq. (23) including its free term are 
significant. 

After transforming Eq. (23) a design guideline was achieved 
for determining the minimum freeboard height, FB, in assumed 
operational conditions and for an arbitrary value of the hourly 
occurrence rate of water ingress into the holds, nz: 


FB = 20 —%2 -La —a,:CBF—a,:CBA -nz (24) 
Oy 
SUMMARY 


In this paper the influence was presented of ship hull form 
parameters on water ingress into the holds, that made it possible 
the formula for elaborating the minimum value of freeboard 
height for open-top containership. The presented method may 
be used in the preliminary designing of ships of this type. The 
investigations were performed by applying a method complying 
with the concept of ,,goal-oriented” approach to formulation 
of design criteria, proposed by IMO. Such approach made it 
possible to elaborate relatively accurate approximations and 
practical guidelines for designing the open-top containerships 
and provides comprehensive information on merits of 
a considered ship already in the preliminary designing phase. 

The discussed investigations consitute only a part of research 
on influence of ship design parameters on the phenomenon 
of water ingress into the holds. Such investigations should 
be extended to research on influence of ship dimensions on 
water ingress into the holds. It may be expected that after 
determination of crucial ship dimension parameters affecting 
the water ingress into the holds it would be possible to make the 
coefficients in Eqs. (23) and (24) dependent on the parameters. 
Results of such research may be very useful in solving the 
design problems of open-top containerships. 

As results from [1, 2, 3, 4] the linear relation between the 
occurrence rate of water ingress into the holds, nz, and the 
volumetric rate of water ingress into the holds, Vn may be 
applicable: 

nz=0: V (25) 
gw 
where: 
nz — hourly occurrence rate of water ingress into the holds 
a  — a constant 
Va hourly volumetric rate of water ingress into the holds. 


For determination of a — value to carry out additional 
investigations is necessary. Knowing the a — value, one is 
able, on the basis of Eqs. (23) and (24), to elaborate a design 
criterion for determining the minimum freeboard height, 
which takes into account the phenomenon of water ingress 
into the holds. It is possible that such research may be based 
on determination of influence of exceedance of the freeboard 
height HO [acc. Eq. (1)] on the volumetric rate of water ingress 
into the holds, Vie 


As resulted from the performed investigations, the 
phenomenon of green water ingress into the holds is affected 
by the initial transverse metacentric height. And, the parameter, 
being operational one, depends on ship load conditions and 
hence it can be hardly considered a design parameter. Fig. 7 
presents the relations which take into account influence of the 
metacentric height and freeboard height on water ingress into 
the holds for a ship of constant dimensions. Such relations may 
be used to plan, in the ship operation phase, a ship loading state 
in advance of ship departure. 


nz [1/h] 


1:25.15 2 


0.3 04 05 0.55 06 0.65 0.7 0.75 08 09 1 


GM [m] 
00-25 25-50 50-75 75-100 W 100-125 M 125-150 


E 150-1750) 175-200 W 200-225 M 225-250 E 250-275 
Fig. 7. Statistical hourly occurrence rate of water ingress into the holds, nz, 
in function of the initial transverse metacentric height GM 
and the freeboard height FB; for open-top containership of the length 
L = 144 m, at speed v = 10 m/s, wave encounter angle p = 60°, 
and significant wave height H, = 8.5 m 
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Finite element local analysis of wave slamming 
on offshore structure 


Bartłomiej Żyliński, M.Sc. 
West Pomeranian University of Technology in Szczecin 


ABSTRACT 


Offshore platforms are exposed to waves slamming event. Waves hitting the columns with 
a high velocity are in many cases the design criteria for column structure. This paper 
focuses on the analysis of wave slamming on floating platform column. Significant for 
wave slamming pressure is load history, which is usually based on model test. Wave 
slamming loads were defined on all four walls of column to assess the worst place. For 
south wall of column three positions exposed to slamming loads between elevation 21.000 
(SWL) and elevation 35.500 were checked. Dynamic analysis has been performed with 


nonlinear FEM program ABAQUS /explicit. The steel was modeled as an elastic-plastic material with 
isotropic hardening. 


Key words: explicit analysis, FEM, offshore structures, wave slamming 


INTRODUCTION 


Wave forces on offshore structures have very often been 
analyzed especially for slender cylinder and walls. For non — 
breaking waves on vertical cylinders those forces are described 
by well known Morrison equation [1]. Morrison equations are 
widely used in engineering practice although can not be used 
for all wave conditions. 

Forces responsible for breaking waves, (especially plunging 
waves) can be more than two times higher than those of non- 
breaking waves of comparable size [2]. Duration of these 
impact forces (also called slamming) is extremely short. Many 
laboratory and experimental studies [3] for regular and non 
regular waves describe intensity, peak pressure, time history 
and distribution on cylindrical surface of piles. Theoretical 
description of impact (forces) was first developed by von 
Karman [4]. 

According to these theories the maximum line forces on 
vertical piles can be given by following equation: 


f =C p Rv’ (1) 

where: 

C, — slamming coefficient, is equal m (due to von Karman 
theory) or 2m (according to Wagner), 

p — mass density of fluid, 

R — radius of the cylinder, 

v — velocity of the mass of water. 


Time history of impact loads was usually calculated 
according to Goda [5]. The new approach to description of load 
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time of impact forces was developed by Wienke [6]. Numerical 
simulation of slamming based on CFD and comparison with 
results obtained from [6] can be found in [7]. 

Complexity of real wave breaking and different types of 
plunging breakers (i.e. early, perfect breaking) maximum forces 
calculated according to theories and measured in experiments 
are well described in [8]. 

Due to this complexity maximum impact pressure on walls 
(as on cylindrical structures) forces calculated according to 
Karman theory (water - hammer analysis) are 8-10 times 
greater than the measurements. Also, other theoretical attempts 
for prediction of maximum impact pressure give unreasonable 
results [8]. 

The safe and economic design of offshore structures such as 
(floating platforms) depends significantly on the designed wave 
loads. In DNV recommended practice [9] prediction of wave 
impact on plates is based on water-hammer analysis models 
but experiments are also recommended, in order to give correct 
estimates of impact loads. 

Space averaged slamming pressure over broader area should 
be calculated from formula: 


(2) 


Slamming coefficient for flat panels should not be taken as 
less than C,, equal 27. There is not given any other information 
about loading time history, area, etc. According to formula (2) 
calculated value of averaged pressure can be very high. For 
ultimate load state (ULS) local plastic deformation of structure 
should be considered. 


] 
P= 5PC 


WAVE SLAMMING ON COLUMN 


Analysed platform has almost identical columns. Wave 
slamming analysis is based on one of columns. The column is 
located at the south / west side of the platform hull. Platform 
north coincides with geographical north for all analysis purposes 
(i.e wave slamming directions). The column structure itself 
shows lack of symmetry. Wave slamming can occur on one of 
four walls. Loads on the inner column walls were reduced by 
20% due to shielding effect from the other columns. 


Sager ot 


seLi0s0 


Elev. 17.500 


L m 


Fig. 1. Model of the analyzed sector of the column 


Analysis methodology / design criteria 


This section describes the methodology used to verify the 
structural strength of column exposed to breaking waves based 
on waves taken from Metocean [10]. Load time history (Fig. 2) 
is derived from model tests and applied in dynamic analysis 
with non-linear FEM program ABAQUS. The slamming event 
has a duration of 160 ms with a peak value at 30 ms. After the 
peak value the pressure decreases to 10% of the maximum 
value. This load history is very similar to the one shown in [8]. 
In the analysis, it is assumed that the main load is the wave 
breaking pressure and platform has no initial speed. 

P 


Prax 


0.3 Pow 


0.1 Prax 


20 ms 10 ms 40 ms 90 ms 


Fig. 2. A design load history for slamming loads on the column face. 


The formula for calculating the pressure is shown below: 


p=LF-+-p-Cyv? (3) 


Slamming coefficient Cs is equal to 27 and load factor LF 
is assumed as 1.30. 


Used slamming coefficient, load factor, wave velocity 
and pressure_are calculated according to DNV recommended 
practice [9]. Due to the mix of air and to the water density p 
is equal to 0.9 t/m?. 


Tab. 1. Wave slamming pressure parameters 


Slamming load on panel (according to DNV CN 30.5) 


100 year RP Hs 
Wave height 
Hb = 1.4 x Hs 
Elevation for top 
of wave crest 
Crest above SWL 
Height of window 
Velocity [m/s] 
Pressure [MPa] 
Line load [N/mm] 


Wave from direction 


> 
S 
= 
~+ 
= 
S 
zZ 


East 90° 


South 180° 


West 270° 


Slamming can only affect areas above still water level 
(SWL), SWL is at elevation of 21 metres. 

Exposed area could be anywhere between still water line 
(el. 21 m) and the top of wave crest. Different positions are 
analysed for south wall of column c1 where calculated wave 
slamming pressure was the highest. Areas used to apply wave 
slamming pressure are shown on Fig.3. Width, height, and wave 
slamming pressure value for each area can be seen in Tab. 2. 
Horizontal width: the flat surface between bilges, but limited 
to a width of about 8 meters corresponding to a sector of 45° 
(with origin in the centre of the column). Vertical height: one 
fourth of the wave height (given in Tab. 2). 


Tab. 2. Areas exposed to wave slamming pressure on south wall of column 


Flat plate pressure 
distribution area Al A2 A3 
south wall of c1: 


Top of exposed area | el. 32250 | el. 29250 | el. 26250 


Width [m] 


Height [m] 


Peak Pressure [MPa] 


Design criteria were defined as follows: stiffeners / web 
frames should have ultimate strength to sustain slamming 
pressure from breaking wave with a return period of 100 
year including a load factor of 1.3, It must be ensured that 
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Fig. 3. Column c1 south wall. Wave slamming pressure exposed area Al -A3 


the structure hit by the breaking wave will stabilize in a load 
history, non-linear analysis using plastic strain formulation 
(acc. to NORSOK N-004 [11]). 


DYNAMIC FEM ANALYSIS 
Explicit solution method 


Abaqus V6.7 explicit was used for the calculations. The 
explicit method is very well suited for the analysis of high 
speed dynamic events and complex and highly non-linear 
problems [12]. 

Explicit calculations are conditionally stable and the 
maximum time increment is limited by the element size and 
the element mass. Automatic time incrementation was used and 
time increment was about 7.75E-06 s. Total time of explicit 
dynamic analysis in ABAQUS was set to 0 - 1.1 [s] and it is 
much longer than the duration of slamming loads (0 - 0.16 [s]). 
The averaged slamming pressure was applied to the structure, 
load history was defined by using amplitude options. Amplitude 
is varying between 0.0 and 1.0, (see Fig. 2). Mass proportional 
damping value was defined as 0.02. For improving numerical 
solution artificial damping is also included in explicit calculation 
procedures by default. Upon wave impact, kinetic energy of 
wave slamming is absorbed by elastic-plastic deformation of 
column structure. 


Material model 


All plates and stringers in the column are made of 420 MPa 
steel. In the calculations, the yield stress of the material was 
taken as 420 MPa. An isotropic hardening model was defined 
in ABAQUS. 

The hardening behavior (engineering as well as true stress- 
strain) of the material is explained in Fig. 3. In a tensile test, the 
force and the elongation are measured. True stress is the force 
per unit of deformed cross section, engineering (nominal) stress 
is the force per unit of the initial cross section. True strain is 
based on the deformed length, whereas engineering strain is the 
elongation divided by the original length. An engineering strain 
of 20% is the same as a true (logarithmic) strain level of 18%. 


Geometry and FEM mesh 


The relevant part of the column for the slamming analysis 
is between the elevation of 17.5 m and 35.5 m. Below deck at 
elevation of 26.5 m, each side of the external shell is supported 
by 2 bulkheads. Over this deck, the shell is supported by one 
bulkhead and one cross-over I-beam from the central void. 
Only the upper part is considered, since the horizontal framing 
is much weaker than below the deck. The outer skin and the 
bulkhead walls are reinforced by stringers HP 320x12 that are 
equally spaced at 0.625 m. In horizontal planes, the outer skin 
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Fig. 4. Hardening behavior of 420 MPa steel 
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Fig. 5. Cross sectional view of the column 


and the bulkhead walls are supported by frames, which are 
equally spaced at 3.0 m. An I-beam connects the corner of the 
central void with the outer skin. 

The FEM model was made with a grid of nodes that were 
spaced in horizontal direction of 312.5 mm and 250 mm in 
vertical direction and the mesh is almost square. In the FEM 
model, the frame flanges were therefore modeled as 300 mm 
high and 20 mm thick. Stiffener web was assumed as 300 mm 
long and 12 mm thick. For the stiffeners on south wall the bulb 
was modeled with T3D2 elements . 

With the given mesh size, there are 5 elements over the 
web of the frames and their behavior will be described quite 
accurately. The stringers have two elements in their web, and 
their behavior is less accurate, but still sufficient. Square shell 
elements of type S4 were used everywhere, except for the 
brackets where some triangular shell elements of type S3R were 
used as well. The brackets between bulkhead and frames were 
modeled as triangles whose short sides are about 1600 mm long. 
The I-beams were modeled with shell elements of type S4. 

Eight millimeters corrosion allowance was used for 
column skin plates (up to 5 meters above SWL) and used 
shell thickness = 16 mm. Thickness of bulkhead and deck 
was 12 mm. 


Load cases 


A wave slamming load is a local load. Waves hit only 
selected (relative small) area of column outer surface. Although 
wave breaking pressures act only local, deformation and stress 
level can be very serious. We can compare this results to results 
obtained from global analysis. Displacements from global model 
(SESAM) have been applied as prescribed displacements in the 
corresponding nodes of the column local model (ABAQUS). 
ABAQUS model used in wave slamming analysis was more 
detailed than the model used in the global analysis. Load cases for 
wave slamming on all four walls correspond to loads described 
in Tab. 1. Load cases on south wall correspond to Tab. 2. Gravity 
was included and gravity constant equal to 9.81 m/s’. 


FEM RESULTS 
Global and slamming loads 


Analysis show that stress level from global model loads is 
very small especially if we compare it to the results from wave 
slamming analysis. Von Misses stress distribution on east and 
north wall of column can be seen on Fig. 6. 

Values of wave breaking pressure acting on south, west, 
and north walls are considerable. Von Misses’s stresses level 
exceeded yield limit of 420 steel. On east wall von Misses’s 
stresses level didn’t exceed yield limit and structure of east 
wall deforms only elastic. For south, west and north walls the 
highest plastic equivalent strains (PEEQ) occurred in the web 
of frame (Fig. 7). The highest strain level was 2.8 % on south 
wall in the web of the frame. Plastic equivalent strains for other 
walls and other parts of column are shown in Tab. 3. 
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Fig. 6. Von Misses stress distribution from global model load [Pa] 
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Fig.7. PEEQ strains on south wall of column 


Tab. 3. Maximum plastic equivalent strains PEEQ in columns parts 


Maximum equivalent plastic strain (%) 


slamming pressure 


(urs) | skin 


south wall 


bracket /frame bulkhead/ 


frame stiffener 


web | flanges | midspan end 


west wall 


north wall 


east wall 
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Plastic strains as a function of time CONCLUSION 


It is important to check whether plastic strains level reminds The results shown herein presents the dynamic nonlinear 
constant after first occurrence. It is possible when dynamic/ elasto-plastic analysis of structure of floating platform. Slamming 
explicit analysis time is much longer then duration of slamming pressures calculated according to DNV recommended practice 
pressure load. Analysis time was set from 0.0 to 1.1 s. Local for ultimate load state (ULS) may compound the difficulties in 
plastic strains were checked for different exposed areas (acc adequately design. Stresses from global analysis are very small 
to Tab. 2) on the south wall. especially if compare to wave slamming results. In the column 
+4.485e-02 structure loaded with the 100 year wave slamming loading 
#4,111e-02 (used load factor 1.3) local yielding occurred. Maximum plastic 
+3, 7776-02 equivalent strains (PEEQ) level was about 2.8%, found at time 
Ei equal to 30 [ms] (maximum of load amplitude). During rest 
+2,616e-02 of analysis time (30 to 110 [ms]) PEEQ strains level remained 
+2,.2426-02 constant. Local plastic deformation of structure should be 
+1,.36%e-02 ; i A 
+1 4956-02 acceptable. The numerical results satisfy equilibrium and 
+1.121e-02 symmetry (Fig. 7). Further investigation of other aspects of 

analysis i.e: local mesh refinement, possible imperfection and 
local buckling of stiffeners can be also quickly performed with 


+7.47S2-03 
#3.737e-03 
ABAQUS system. 


+0.000e+00 
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Failure criteria of viscoelastic materials 
under multiaxial loads 
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ABSTRACT 


The load capacity of homogeneous, isotropic viscoelastic materials subjected to multiaxial 


stresses is considered. For this purpose uniaxial equivalent stresses in selected load cases 
are determined and relevant criteria applied. It is shown that in the case of constant 
load the yield criterion does not differ from that for perfectly elastic materials. Similar 
conclusion has been drawn for the fatigue and yield criteria at in-phase and out-of-phase 
stresses. On the contrary, the criteria derived for viscoelastic materials subjected to 
periodic loads do not coincide with those for perfectly elastic materials. 


Keywords: viscoelastic material, multiaxial stress, fatigue criteria, yield criterion, static load, vibratory load 


INTRODUCTION 


Even in the region below the limit of proportionality, metals 
are not perfectly elastic, and their deformation is accompanied 
by internal friction. As a result, their strain response to static 
loads is time-dependent, and when the load varies in time, 
the strain lags behind the stress. In particular, after the rapid 
shortening of a bar that occurs when its axial tensile load is 
removed, the bar continues to shorten gradually until the initial 
length has been reached. This gradual return to the initial 
length following unloading is called anelastic behaviour or 
viscoelasticity [1, 2]. Anelastic strain differs from plastic strain 
because it is recoverable rather than permanent, and it differs 
from elastic strain because it is recoverable at a rate which is 
slow in comparison to the rate of recovery of elastic strain. The 
elastic behaviour of metals and structural steels is much more 
significant than their anelastic properties, which enables the 
strength criteria to be formulated with the aid of the Hooke’s law 
for perfectly elastic materials. However, as shown in the present 
paper, dissipative properties of engineering materials may be 
also taken into account in design considerations, especially in 
the case of periodic load when an equivalent reduced stress is 
to be determined in calculations of a structural limit state. 

During the last two decades, the emphasis in structural 
design has been moving from the allowable stress design to 
the limit state design [2]. Generally, four types of limit states 
may be specified: 
serviceability limit state, 

x ultimate or yield limit state, 
x fatigue limit state, 
x accidental limit state. 


x 


This paper concerns the yield and fatigue limit states of 
viscoelastic materials under combined static and vibratory 
loads. 


YIELD CRITERION OF VISCOELASTIC 
MATERIALS UNDER MULTIAXIAL 
STATIC LOADS 


The stress and strain at which a material either begins 
to yield or fractures due to a uniaxial load can be measured 
relatively easily. But for an arbitrarily shaped body under 
arbitrary loads, the prediction of yield or fracture is very 
difficult. Some criterion is needed to make predictions without 
testing every material under every possible loading. An ideal 
criterion would be one that is based on a simple uniaxial test. 
Then the normal stress, normal strain, shear stress, the strain 
energy, or the distortion energy, among other possibilities, 
could be taken into account. Each of these reaches its failure 
value at the same load in a uniaxial test, but this is no longer 
true if the state of stress is either two- or three-dimensional. 
Therefore various theories regarding the initiation of yielding 
have been developed. The earliest, and the simplest, relation 
describing the conditions for initiating plastic flow under static 
or quasi-static load is the shear stress law. This law states 
that the metal will yield when the largest shear stress reaches 
a critical value, irrespective of the stress state. A somewhat 
more accurate law is the “energy-of-distortion” criterion, also 
called the Huber-Mises-Hencky (HMH) strength theory [2-4]. 
It embodies the physical hypothesis that yielding occurs when 
a certain critical value of distortion energy is reached, which 
corresponds to that where the equivalent reduced stress, © 
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reaches the yield strength, R,, determined from the uniaxial 
tension test. Experiments show that this is an excellent criterion 
for the yielding of ductile materials [2]. This criterion is also 
useful for ductile materials under multiaxial proportional 
loading and high-cycle fatigue [5]. Therefore it will be also 
used in the present paper. 

To begin with, let us consider a perfectly elastic solid under 
general state of static load resulting in normal and shear stress and 
strain components ©.» £. Tg» Yk G = X, Ys Z; K = xy, yZ, ZX). 


jo? “jo” 


The elastic strain energy per unit volume is [4]: 


1 
WV, = Zoe ‘Tete (1) 
j k 


where: 
by, = =F -vlo,, +0, )| 
1 
e, = =b, vlon +o, )| P 
1 
Eo = aA -vlo +o, )| 
Yio = G 


In Eqs (2), expressing the Hooke’s law for multiaxial stress 
in elastic solids, E is the Young modulus, v is the Poisson’s 
ratio, and: 

E 


j 2(1+v) G3) 


is the shear modulus. Substitution of Eqs (2) and (3) into Eq. 
(1) results in: 


1 > 
Fu +O, +6,,) + 
y 


3 + Tyzo + Tao ~ O oO yo + 


(4) 
+(1+v) ii 


ol 
o7 E 
E G9 2 = O w0 xo 

which can be rewritten as [4]: 


Y, = Yao T Wo (5) 


lo. a i + 


_l+v ( y 
Wao ~~ 6E + Oso -0 


where: 


+(6,,-O) +] © 
+ 6(t2,, +, + t2.) 


is the strain energy of distortion, and: 


Oo +0. tO, ) (7) 


is the strain energy of volume change. According to the 
HMH theory, at yield, the distortion energy at a point in a three- 
dimensional state of stress is equal to the distortion energy at 
yield in the uniaxial case. Thus: 


Wiss = Yao (8) 
where: 
Y = lto o? (9) 
deo 3E eo 


14 POLISH MARITIME RESEARCH, No 1/2009 


is the strain energy of distortion per unit volume under reduced 
stress o „İn tension. Hence [2, 4, 5]: 


o, +07, +0, =O oyó =O 
Oo = > 5 5 (10) 
z 9,9 x5 + 3(c2,, + Tyo + 2.) 
and the criterion in question reads: 
Oa I R, (11) 


Now suppose that the considered load is applied to 
a viscoelastic solid at the time t = 0 and the stress components 
(o and t,, are maintained constant. Then the constitutive 
equations for strains are [6]: 


E, 
e,(t)= E| l-e " 


i (12) 
y,(t)= ol Ie” 
and the elastic strain energy becomes time-dependent: 
w(t)= > ze jo&(t)+ 2 To¥e(t)] 03) 


In Eqs (12), n is the coefficient of viscous damping of 
normal strain in the Kelvin-Voigt’s model of the material 
[2, 6]. Substitution of Eqs (2) and (12) into Eq. (13) gives: 


E 
-=t 


w(t)=w,| 1-e” (14) 
with: 
lim w(t) =W, 


to 


(15) 


Hence it is clear that Eq. (10) and the criterion (11) are 
applicable also to viscoelastic materials under constant loads. 


FAILURE CRITERIA OF VISCOELASTIC 
MATERIALS UNDER MULTIAXIAL 
HARMONIC LOADS 


Under dynamic loading conditions, the strength of 
a structural element is degraded due to the cyclic application 
of load or strain which may lead to fatigue damage. It is why 
in this paper not only the yield criterion but also fatigue criteria 
are considered. 

The fatigue load that a structure can withstand is often 
significantly less than the load which it would be capable of if the 
load were applied only once. For the case of combined dynamic 
loading, where torsion, bending and/or tension loads vary in time, 
no one strength theory is universally accepted, and all existing 
multiaxial fatigue criteria can demonstrate large scatter [7]. 
Only multiaxial in-phase stresses with constant principal 
directions can be treated fairly well using the conventional 
strength theories [8]. In [9] an attempt was made to extend the 
application range of the HMH theory to other stress states, in 
particular to those with variable principal directions. In what 
follows, the results obtained in [6, 9] are taken into account. 

If a multiaxial harmonic load is producing in a perfectly 
elastic solid the zero mean stress with in-phase components: 


o. =0. snot: j=X,y,Z 
j ja J y (16) 


Tk =T Sinot ; k=xy,yz,zx 


the equivalent reduced stress in tension-compression takes the form: 


o, =6,, sinwt (17) 
where: 
Oja Tka — amplitudes of the stress components ; œ — circular frequency ; o „— amplitude of the reduced stress given by [8]: 
? 2 2 1/2 
Ori +O), tO =O, Saa" S o a a 
Oa = i noe (18) 
=0..6. _ +3(¢2 cok +T +T) 
Consequently, the criterion of an infinite fatigue life and the yield criterion read: 
Oe s Zi (19) 
Oa L R, (20) 
where: 


Z „~ the fatigue limit of the material under fully reversed tension-compression. 


In the high-cycle fatigue regime, the criterion of a finite fatigue life based on the S-N curve (Wöhler’s curve) is commonly 
accepted [10]. According to this approach, in the case of fully reversed tension-compression under the stress (17) one gets: 


Na om 
—o,, <l (21) 
K 
where: 
N, — required number of stress cycles to achieve a given design life 

K — fatigue strength coefficient 

m — fatigue strength exponent. 
If, however, the in-phase stress A is applied to a viscoelastic material, its strain response is [6]: 


Tba- vlo., +0, „Jbinot-a) 


=k. -v(6,, +O, )|sin(ot-«) 
(22) 


“ret vlo., +O, _)|sin( lot- a) 


_2(l+v) _ + v) 


JE +o? 


So, the elastic strain energy per unit volume: 


l 
-3[Eee Ean] (23) 
j k 


(Q) 
, sin (ot - a), a= aretg— 


becomes: 
y= ———— ,, lo. - ul, TO- Jlo, [o,, -v(o,, +6, )}+,, Is, -vlo +0, )|+ 
= +n @ 
>|. ‘ ] 1 3 (24) 
+ 2(1+v)>° ta psinotsin(@ot-a)= ibe, +6,,+6,, P + 
k JE? +n°o? L2 
+(1+v)(c2,, +1 Tja +T, =O aO ya TO paO na -0,5 ) sino tsin(o t-a) 
On the other hand, for the ne stress (17) and corresponding reduced strain: 

E, po a sinlot-a) 025) 

the elastic strain energy per unit volume is given by: 
W. = 0.8, “Sr sing tsin(at — a) (26) 
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Consequently, the relationships for distortion energies per unit volume 
in the actual and reduced stress states read: 


Y= ze -Oy, P +o, -0n i +(c,-c,,) +6(t2,, +T HTa )sino tsin(o t—a) (27) 
6E? +07 


l+v 2 . 
s =— o sinotsin(ot-a) 
3E? +70? 
Equating the right-hand sides of Eqs (27) and (28) yields again Eq. (18) 
so that the criteria (19) through (21) remain valid. 


In the case of non-zero mean stress with in-phase components: 
o. =0. +0. sinoat ; tT =T, +T, sinwt 
J Jo Ja k ko ka 
the equivalent reduced stress can be expressed as: 
o, =0_ +o_ sinat 
e eo ea 
where: o „is the mean value of the reduced stress. In Eq. (30), there are two unknown parameters 
of the reduced stress (o,, and o) which cannot be determined from the single equation, 
corresponding to the HMH theory, without additional assumptions. On the other hand, 
since this theory can be applied to the stress (29) when: 
Om #0 > 1, #9 ; Oa = Ta = 0 
and when: 
Om = Tim 7 O ; o, #0 > 1,, #9 
it should be also applicable to the stress (29) when: 
o #0 ; 1,40 G o, #0 > Ta #0 


jm 


Therefore the following hypothesis (“average-distortion-energy strength hypothesis”) was formulated [11]. 


The reduced stress (30) is equivalent in terms of static and dynamic effort 
of a material to the stress (29) if: 
(i) the time-independent parts of distortion energies per unit volume 
in these both stress states are equal 
(ii) the reduced stress and the stress components (29) 
have the same frequency 
(iii) the integral time averages of instantaneous values of distortion energies per unit 
volume in these both stress states are equal. 


The distortion energy per unit volume in the general state of stress is given by [4]: 


Wa 40.6. ~¢,,)+6,(¢, ~¢,,)+6,(€, -2,)+ Dun | 
k 


where: 
l 
En =5. +e, +e,) 


When the stress (29) is applied to a viscoelastic material, according to Eqs (2), (12) and (22) one gets: 


E, = Ii, -vlo +0, )] l-e" — h, -vlo., +0, )sin(ot-a) 
E i E’ +170" ' 
E 
E, =—|c,, vlo, On )| e" |b -v(o,, +6,, )|sin(ot —c) 
E E’ +o? 


o 
+- tm sinfot-a), a= arctg T 
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(28) 


(29) 


(30) 


(31) 


(32) 


(33) 


or, after sufficiently long time: 


e =b,- vlo, +0, Jess 


e, = Ii., -vlo +o,,) +o, -vlo +6,,)}sin(ot-«) 
yY Re E? 22 6) 


e.- sk. -0 i 


ane 


lo.. - vlo., +0, )[sin(ot —a) 


lo.. -vlo,, +6,,) sinlot-a) 


1 
Oo aÅÃ—— 
i 
t,, sin(ot—«), a = arctg 22 
ae +o? E 


Substitution of Eqs (3), (32) and (34) into Eq. (31) results in: 


Gao 
y= Ozo Oyo Om -O x00 yo -O p00 — Fn 0 +32 Tat 
3E 
l+v sÀ 9? 2 2 . . 
+———_—| 62, +0, +07, - 6,0), —5,,5,, —5,,0,, +3>_%, |sinotsin(at—a)+ 
34E? +o" l l l k 


hbo, }+0,,(2 }+0,,(2 ) 
+ SE Oya ZO xo Oyo -S Oy, Oo O45 —O8,,/17 95, O — Oio — Oyo +F 


l+v 


+6% tuti sin ot+— lb. 2o. —6,, ~o,,)+ olo, -0a -0,)+ 
k 6E? +07 
5, (26,, -0 TAn 6> tutu [infot -a) 
k 


Hence the distortion energy per unit volume of the viscoelastic solid under the reduced stress (30) is: 


+0 > Il+v 
+ 


l 
W e = 5 Gaos 4 
“  3E 3E? +o? 


l+v l+v 
+——6,,6,, sinoat + — 0 0 sin(ot — a) 


3E HWE tno a À 
The condition (i) gives Eq. (10), and the condition (iii), i.e., 


E t 27 
= fvadt = = give ; T=— 
a) 


o? sinot(@t—a)+ 


T 0 
leads to Eq. (18). Then the criterion of an infinite fatigue life reads [2, 10]: 
ie) O a 
eo + ea < 1 
R, re 


and the yield criterion is: 


m 
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(35) 


(36) 


(37) 


(38) 


(39) 


(40) 
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Eq. (40) is valid if: 
(0) 
Z. <——— <L 


j (41) 
l aps eo 
where: 
L- the maximum stress amplitude under fully reversed tension-compression 
in the high-cycle fatigue regime (above which the low-cycle fatigue may occur). 
Another stress state which can be dealed with the aid of the aforementioned hypothesis 
is that with non-zero mean out-of-phase components: 
5,=6,,+6,,sin(at+@) ; T= To Ta sin(at + Q) (42) 


where: ọ, and ọ, are the phase angles. With the known strain response of viscoelastic materials 
to the stress (42) [6] and with the reduced stress in the form (30), it is easy to prove that the mean value 
of the reduced stress is given again by Eq. (10), but now the formula for its amplitude reads: 


5 > 2 1/2 
(o Ta i o a a 0 a 2 A o coslo, -~ 0, )- 6,5, coslo, -ọ,)+ (43) 
Oa = 


ea ? R ? 
97,9 x4 cos(¢, E Px )+ C Caa + Tyza + Taxa ) 


It means that in the considered above load cases the fatigue and yield criteria for elastic and viscoelastic solids coincide. 
Of course, to be on the safe side, in Eqs (38) through (40) arbitrary safety margins may be introduced. 


FAILURE CRITERIA OF VISCOELASTIC MATERIALS 
UNDER MULTIAXIAL PERIODIC LOADS 


In general, structural elements and machinery details are simultaneously subjected to static and dynamic loads. In addition 
to the time-varying stress, the steady stress resulting from static load, and the mean stress (the average of the maximum and 
minimum of the cyclic stress) influence the strength of an element. It is important to recognize that the total strength of the 
element is altered if residual stresses (caused by cold forming, heat treatment, welding, etc.) exist. Since residual stresses have 
a similar influence on the fatigue behaviour of materials as do mechanically imposed constant stresses of the same magnitude 
[12], in what follows no distinction will be made between any kind of static stresses. 

Among the types of dynamic loads encountered in practice, one of the most important is periodic load. In the general state 
of periodic stress, its components can be expanded in Fourier series: 


6, =6,,+)> On sin(not +9) ; j=xX,y,Z 
g (44) 
Ty = T +È Tu Sin(not+ Pn) ; k= xy, yz, zx 
n 
where: 


O; Tko — mean values ; @,,, Q,,, — phase angles of n-th terms ; œ = 27/T — fundamental circular frequency ; T — stress period. 


jo” 
Our aim is to determine such a reduced stress: 
o, =O + 6,, Sinot) (45) 


of mean value o „amplitude o,, and circular frequency œ „ that would be equivalent to the original stress in terms of fatigue and 
yield strengths of viscoelastic materials. For this purpose the theory of energy transformation systems [13] can be used, which 
links the dissipated energy with the breakdown time of a system. According to this theory, two stress states are equivalent in 
terms of time to failure if the energies dissipated internally and externally in both these states are respectively equal. Apparently, 
the dissipated energy can be estimated by evaluation of certain symptoms, e.g., the externally dissipated energy by the vibration, 
and the internally dissipated energy by the temperature [13]. For viscoelastic materials the use can be made of the following 
relationship derived in [14] for the internally dissipated energy per unit volume under periodic stress with the components (44) 
during the period T: 


nsina, 


l 5 i 5 
o(T)= oT} ——— 6, +05, +07, ~2v|c,,,, coslo., =op) 
2 5 VE’ +(nno)’ (46) 


+O nOn cosl yn Pm )+ Opon COS(P,, — Pye )|+ 2(1+ vl, HT + Then } 


where: 


no 
t= arct 7° (47) 
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Consequently, the following equivalency condition can be postulated: 


p(T) = (T) (48) 
where: 
¢(T) — the energy dissipated in the viscoelastic material per unit volume 
under the reduced stress (45) during T seconds. Under assumption that: 


oO, =ro (49) 
where: 
r—natural number to be determined, on the basis of Eqs (46) and (47) one can write: 
l sina, A 
o.(T)= Pa T+, (50) 
E’ +(nro) 
ro 
Q, = art e (51) 
E 
; Thus 
rang > 
-o =A 
E’ +(nro) (52) 
where 


A= yia p +0, +07, ~2v]o,,6,, coslo., E 
» E? +(nno) (53) 


+6,,0., coslo., -Pn )+ OnO COSlOn — ,, )|+ 2(1+ v)(t?, +t HTa } 


As for the remaining equations, necessary for determination of the reduced stress (45), in [9] it was assumed that the energy 
dissipated externally by materials subjected to multiaxial loads is proportional to the average strain energy of distortion per unit 
volume. Bearing in mind the role of distortion energy in evaluation of effort of ductile materials under multiaxial stresses, in 
the present paper the same assumption is retained. Derivation of similar equations based on the strain energy of volume change 
or total strain energy is analogous. 

In conformity with the aforesaid, we have: 


1? 1? 2T 
—|w,dt=—lw, , T=— (54) 
7! Vide al Pa i 
where: 
+0 5 Il+v 2 : 
Y =—_9., +O sinro tsin(re t — a )+ 
de 3E eo > 2 ea e 
3E’ +(nro) 
l+v l+v (53) 


+—— O aO SINTO t + G,,,6,, Sinfrot-a,) 


3E 34E? +(nroy 


is the strain energy of distortion per unit volume under reduced stress (45) and y, is that under periodic stress (44). To determine 
the latter from Eqs (31) and (32), the strain response of a viscoelastic material to periodic stress must be known. This problem 
was solved in [6] to give: 


e=e,+> H, Imla, el net an | (56) 


where: 


T 
£, = le. Eyo Ea Yxyo Yyzo Yzxo | 


is the vector of mean strain components given by Eqs (2), 


1 -v -v 0 0 0 
-v 1 -v 0 0 0 
1 -y -v | 0 0 0 
H, B 
E2+(nnoy| 9 O0 0 2Al+v) 0 0 
0 0 0 2Al+v) 0 
0 0 0 0 0 2(1+v) 
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is the matrix of dynamical flexibility of the viscoelastic material at the load circular frequency no, 


Gy = [BFF Ton A | 


xyn “yz * zxn 
is the vector of complex amplitudes of n-th terms of the stress components, defined as: 
= — Qn. F _ iPkn 
On = Fin > Tin = ne " 
iis the imaginary unity, Im is the imaginary part, and: 
nno 
a, = arctg— 
E 
is the phase angle of n-th terms of the strain components. 
Eqs (31), (32) and (54) through (56) lead to: 


I+v , (1+ v)cosa, 2 l+v 
> Ox m = On = 4 

3E 64E? +(nro)’ 3E 
+T, HTI ae ania ee eaten (IE Š +0, +0, —O nI yn coslo., -9 )+ (57) 


n yE? +(nno) 


2,2 2 2 
=O nOn coslo., Haa coslo, -Pa )+ 3(r?,, PTa T Tan )] 


Hence: 


z 
[o?, +05, tO, Om TT t 


=o. (+3(t? 


xyo 


a. =O: o + D5, +07 o o a S 40. SE e E s +3h? S tea a T, a (58) 
which results in Eq. (10) for the reduced mean stress, and: 


cosa, 
E +(nro)’ 


where: 


Oa =B (59) 


B= y-o +65, tOn O mO yn coslo., ~ yn) Tyg Fy cosl mPa )+ 
n yE? +(nno) (60) 


Ty n 


From Eqs (52) and (59) one obtains: 
1/4 


= H eliro) (42 +1°B’) (61) 
r? 
and: 
rtga, = = (62) 
° B 
that is: 
‘NO A 
r“ = =— (63) 
E B 
Hence the real number p , close to the natural number r, 
r = Round(p) (64) 
is: 
= -= É (65) 
j noB 


With these results, the criterion of an infinite fatigue life and the yield criterion can be expressed by Eqs (38) and (39), 
whereas the criterion of a finite fatigue life becomes: 
m 


ca 


roT,| o <1 
2nK Oo (66) 


where: T , — the required design life. 
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CONCLUSIONS 


O On the basis of the Huber-Mises-Hencky (HMH) theory 
and constitutive equations for strains in viscoelastic solids 
subjected to multiaxial static loads, it is shown that the yield 
criterion does not differ from that for perfectly elastic solids. 

O On the basis of the HMH theory and constitutive equations 
for strains in viscoelastic solids subjected to multiaxial zero 
mean in-phase loads, it is shown that the criterion of an 
infinite fatigue life, the criterion of a finite fatigue life and 
the yield criterion do not differ from those for perfectly 
elastic solids. 

O With the aid of the average-distortion-energy strength 
hypothesis and constitutive equations for strains in 
viscoelastic solids subjected to non-zero mean in-phase 
loads, as well as subjected to non-zero mean out-of-phase 
loads, it is shown that the aforementioned criteria do not 
differ from those for perfectly elastic solids. 

O On the basis of the theory of energy transformation systems 
and constitutive equations for strains in viscoelastic 
solids subjected to periodic loads, it is shown that the 
aforementioned criteria do not coincide with those for 
perfectly elastic solids. 

O The load capacity of a homogeneous, isotropic viscoelastic 
material at a given temperature is completely defined by 
the Young modulus, tensile yield strength, S-N curve for 
tension-compression, Poisson's ratio and coefficient of 
viscous damping of normal strain. 


NOMENCLATURE 

A,B — quantities defined by Eqs (53) and (60) 

E  — Young modulus 

G  — shear modulus 

H, — matrix of dynamical flexibility of the viscoelastic material 
at the load circular frequency no 

i — imaginary unity 

Im — imaginary part 

K  — fatigue strength coefficient in equation of the S-N curve 
for tension-compression 

L  — maximum stress amplitude under fully reversed tension- 


compression in the high-cycle fatigue regime (above 
which the low-cycle fatigue may occur) 


m  -— fatigue strength exponent in equation of the S-N curve for 
tension-compression 

n — natural number 

N, — required number of stress cycles to achieve a given design life 

r — natural number given by Eq. (64) 

R, — tensile yield strength 

t — time 

T  — stress period 

T, — required design life 

Z, — fatigue limit under fully reversed tension-compression 

a — phase angle of the strain components under in-phase loads 

a, — phase angle defined by Eq. (51) 

a, — phase angle defined by Eq. (47) 

Y, — k-th shear strain component (k = xy, yz, zx) 

Yo ~ k-th strain component under static load, mean value of 


k-th strain component 
€ — vector of the strain components 
€ — reduced strain 
E — j-th normal strain component (j = x, y, Z) 
€ — j-th strain component under static load, mean value of j-th 
strain component 
quantity defined by Eq. (32) 
vector of the mean values of strain components 
— coefficient of viscous damping of normal strain 
— Poisson's ratio, 
number given by Eq. (65) 
— reduced stress 
— amplitude of the reduced stress 
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— reduced stress under static load, mean value of the reduced 
stress 

— j-th stress component 

— amplitude of j-th stress component 

— amplitude of n-th term in Fourier expansion of 9, 


o- — complex amplitude of n-th term in Fourier expansion of 9, 


— j-th stress component under static load, mean value of j-th 
stress component 


6 — vector of complex amplitudes of n-th terms in Fourier 


expansion of the stress components 

— k-th stress component 

— amplitude of k-th stress component 

— k-th stress component under static load, mean value of 
k-th stress component 

— amplitude of n-th term in Fourier expansion of t, 

— complex amplitude of n-th term in Fourier expansion of t, 

— phase angle of j-th stress component 

— phase angle of n-th term in Fourier expansion of 9, 

— phase angle of k-th stress component 

— phase angle of n-th term in Fourier expansion of t, 

— dissipation energy per unit volume 

— dissipation energy per unit volume under the reduced stress 

— strain energy per unit volume 

— strain energy of distortion per unit volume 

— strain energy of distortion per unit volume under the 
reduced stress 

— strain energy of distortion per unit volume under the 
reduced static stress 

— strain energy of distortion per unit volume of perfectly 
elastic solid under static load 

— strain energy per unit volume under the reduced stress 

— strain energy per unit volume of perfectly elastic solid 
under static load 

— strain energy of volume change per unit volume of 
perfectly elastic solid under static load 

— circular frequency, fundamental circular frequency 

— circular frequency of the reduced stress at periodic loads 
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Hydrodynamic performance 
of coaxial contra-rotating propeller 
(CCRP) for large ships 
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ABSTRACT 


This paper describes a coaxial contra-rotating propeller (CCRP) system to calculate the 
hydrodynamic characteristics and then obtain the optimum operational condition to install 
on two different large bulk carrier and VLCC. The method is based on boundary element 
method (BEM) to obtain the hydrodynamic performance of any complicated configuration 
such as CRP system, and then the optimum propeller data is obtained by the systematical 
method at the design condition. We prepared a software package code, namely SPD 
(Ship Propeller Design), which has model mesh generation, solver and numerical output 


results. The comparison of the propulsive performance was made between the propeller alone and CRP 
arrangement. Major finding include good agreement between predictions using the numerical code and 
experimental data for both ships. 


Keywords: boundary element method, contra-rotating propeller, CRP, optimum efficiency 


INTRODUCTION 


The main aims of the marine propulsor designers are to 
increase thrust, diminish torque, improve efficiency and to 
save energy. Most conventional propulsors are working behind 
the ship hull where the flow is non-uniform, unsteady and also 
limitation of the propeller diameter due to ship stern. Single 
propeller produces cavitation at the heavy load condition, 
but CCRP may provide moderate load on each blade and 
prevent cavitation. Another case is that the single propeller 
generates torque, while CCRP cancel the torque and improve 
the efficiency [1, 2]. 

The history of CRP goes back when a patent was applied 
by Ericsson (the inventor in 1836) to 45 feet ship. In 1909 and 
1939, Italian Navy and US Navy had experimented CRP on 
a 46 feet and 70 feet steam ship, respectively. Rutundi [3] made 
a comparison test between CRP and conventional propeller 
for a 3500 tons naval training ship and has reported an 18% 
improvement in the propulsive performance. 

Since then, CRP has well been used for torpedoes, small 
vessels, and of course for aircraft, but there is a difficulty in 
producing a reliable CR shafting which can support the large 
power for application to large merchant ships. In 1988, MHI 
(Japan) has succeeded in retrofitting 4200 GT with a CRP [4], 
and in the same year, [HI (Japan) has completed the shop test 
of Juno’s CRP system at the outset. Having been equipped 
with this CRP system, Juno dealt with the official trial in witch 
she has achieved a 15% power saving. After that, some other 
experimental research have been done at NMRC (in Japan) 
(before name was SRI) by Ukon, [5, 6, 11, 13, 14]. 
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From the numerical approach during two decades, some 
work has been done to obtain a better understanding of 
the system in order to predict its effect on hydrodynamic 
performance. Yang C-J [1] and Islam et al. [8] was used lifting 
surface theory and vortex lattice method to calculate the blade 
loading. 

Recently, more attention is being drawn to the development 
of the contra-rotating podded propulsor (CRPP) system for 
ship propulsion because of its attractively high energy saving 
rate as well as lower cavitation and better hydrodynamic 
performance. 

In the current arrangement, a CRPP is placed at the forward 
end of a pod which is aligned with the local inflow. The 
powering and cavitation experiments show the performance 
prediction agree well with measurement. 

Nishiyama and Sakamoto [10] designed a CRP system at 
IHI and on bulk carrier and VLCC. 

This paper deals with the following subjects: 
numerical Method of BEM 

mesh Generation 

design concept of CRP on Behind of Ship 
design concept of CRPP 

application of the CRP system on the Ship 
systematic design 
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This paper is firstly calculated the hydrodynamic 
performance of the CCRP using BEM. The calculated results 
are predicted well with experimental measurements. Systematic 
design is finally employed to the two ships types (Bulk carrier 
and VLCC) to obtain the design points. 


DESIGN CONCEPT OF CRP 
Procedure of CCRP designing 


With the principal particulars of the hull and the main 
engine given, a CRP can be designed as for any single propeller, 
namely, through reiteration of the following two steps: 

(1) Propeller Design in uniform flow: try to find for an 
optimum design that satisfy the given engine power and 
revolution speed; then 

(2) Modification by considering the non-uniform flow: the 
propeller designed thus is further modified, so as to better adapt 
to the non-uniform flow over what has already been considered 
macroscopically in step (1), in such terms as propeller 
cavitation, propeller-induced vibration, and strength of blade 
for the propeller performance in non-uniform flows. 


Designing CCRP in uniform flow 


There are two approaches in the design of CRP: one based 
on lifting surface theory [9], and other utilizing the design 
diagram based on open water tests on systematic CRP model 
series [5]. More precise results could be expected from the 
systematic series data than from the theoretical approach. 
However, a review of past reports for CRP, the method based 
on series test data can be found. 

Fig. 1 depicts the design procedure and algorithm. Here, 
the relation between ship speed and resistance (R, — V,) for the 
hull, self propulsion test or some empirical formulae (t, w, np), 
the number of blade (Z) and axial distance between two 
propellers (x,) are given. According to systematical method 
and flowchart, first estimate the revolution number of propeller 
and its diameter. Then using numerical SPD code, the 
hydrodynamic characteristics (open water) of propeller alone 
and CRP system are calculated. 

How to calculate the optimum efficiency from the open 
water diagram? This is the systematic method to obtain 
the optimum design of the propeller for the ship. From the 
resistance and the self propulsion test or some empirical 
formulae, the following relation may be expressed: 

* For single conventional propeller: 


= ——; = =] = 
7 pn7D* pD’v, 


(1) 
E -= AJ? 
(1—t)pD* [yva -w 
* For CRP system: 
_ erly — 
pnp Dy 
= —— S = BJ? 
(l ~ t)pD;[V. (1 a w)] 


This K, is quadratic function of J and intersects with open 
water curves of K, — J and the optimum solution is obtained 
from the intersection point, so all the coefficients are obtained 
from this point (Kp K,, n, J). Then, the thrust and torque are 
obtained from the following: 


T=KQ,,,)pwD*;.T 


> 
? ~ Calculated(i+1) — 


Ky 
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Q=K,(,,)pn°D® ; Q,,, <Q 


P 
Ww 
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, _ npQpt+n,Q, 
Ko = 3+ (4) 
pnr D; 
n -J Kr . j=- Va (5) 
= 2nKg üi 


The reiteration is continued unless the thrust should be 
bigger than the required one and torque should be less than 
the previous stage. 


NUMERICAL APPROACH 
Potential based boundary element method 


Suppose the forward propeller of a CRP with Z, blades 
rotates in the left hand (counter-clockwise) direction at 
a constant revolution number n,, while the aft propeller Z, 
blades in the right hand (counter-clockwise) direction at n,, 
and the CRP as a whole advances at a constant speed V,. 

Assuming inviscid, incompressible and irrotational flow 
in the volume around and inside the body, a potential function 
exists for the perturbation velocity b created by the propeller 
movement in the volume which satisfies the Laplace’s equation. 
By applying Green’s theorem for perturbation velocity 
potential @ at any field point on the body surface, we can get 
the following integral equation on the propeller and its trailing 
vortex wake. 


af 1 av(q)( 1 
2nb(P)= Í 46(q)<—|—— |- eak dS+ 
Spl  PaARO;q) Ong \R(p:q) 

ô l (6) 
+ | AAD —— 
Sw Ong \ R(p;q) 


R(p;q) is the distance from the field point p to the singularity 
point q. This equation may be regarded as a representation of 
the velocity potential in terms of a normal dipole distribution 
of strength (p) on the body surface S,,, a source distribution 
of strength 0o/on on S,, and a normal dipole distribution of 
strength Ad(q) on the trailing wake surface S 


Boundary conditions 
The strength of the source distribution in equation (6) 


is known from kinematic boundary condition (KBC) as 
follows: 


-y i=-M +5} (7) 


where: 
n — denotes the outward normal unit vector. 

The strength of dipole distribution is unknown and equal 
to the perturbation potential on the propeller or to the potential 
jump in the trailing vortex wake. On the wake surface S, the 
velocity is considered to be continuous while the potential 
has a jump across the wake. It is expressed in the perturbation 
potential as: 


aa AE 
C 
aap 


ôn Js \@n ôn 
W 


(8) 
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(Ab), =$°-$' =T ©) 


where indexes B and F mean back and face sides of the 
propeller, respectively. 

Another important physical boundary condition is the 
Kutta condition and its implementation. This equal pressure 
Kutta condition is applied to determine the unknown AQ,,, 
of the dipole strength on the wake surface. In the numerical 
calculation, the pressure Kutta condition at the back and face 
surfaces of the trailing edge (TE), can be expressed as: 


(10) 


A direct solution of the resulting system of equations 
obtained from discretized Green’s formula for the perturbation 
velocity potential (6), along with equation (10) is difficult due to 
the nonlinear character of the equation (6) therefore, an iterative 
solutions algorithm is employed to solve the problem. We focus 
on the numerical implementation in the following section. 

Discretization of equation (6) leads to a linear system of 
algebraic equations for the unknown ọ as: 


Ship Resistance at Design Speed: (Rr — Vs) 
Propulsion factors: t, w, Np 


Number of blade: Z 
Axial Distance Between two propellers: Xa 


First Estimation of the Propeller RPM and diameter 
J=V,/nD 


Calculation of equivalent operating 
condition for forward propeller 


Calculation of equivalent operating 
condition for aft propeller 


Convergence of 
equivalent operating 
condition 


Calculation of Optimum propeller efficiency 


Design of forward propeller 


Design of aft propeller 


Based on Potential 
Boundary Element 
Method to calculate 
the CRP Performance 


Intersection of open water diagram 
(K; — J) with K, = AJ- 


A =f(p, Rr, t, w, Vs, D) 


Calculate the Thrust and torque at the 
optimum efficiency 
T=K(J)pn°D* ; K p nD? 


Convergence of Thrust 
T Calculated 2 Trequired 


Fig. 1. Calculation flowchart of CRP 
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(11) 

2 K Ntot Ab 

k | © 
+>] 2D 2s |= 
Ni | k=! j=1 ` (On 
i= 1,2.....Niot ; where: N = 2xNxM 
K - number of blades 

Where Di ; Wi (dipole distributions on body and wake 


surfaces) and oF ' (Source: distribution on body) are influence 
coefficients on panel j j acting on the control point of panel i. 
Those influence coefficients are nearly evaluated analytically. 
The use of quadrilateral surface panels instead of planar panels 
has been found to be important for the convergence of the 
present potential based boundary element method. It is found 
to be especially so when applied to the highly skewed propeller 
and twisted shape. 


Calculation of induced velocity 


From Green’s theorem in the potential field, equation (6), 
we can alternatively construct in the velocity field. Taking the 
gradient of the perturbation velocity potential at any field point, 
the induced velocity which can be expressed as: 


2 é 1 
4nv(p) = + 
P d P ôn, (R(p;q) 
alq) 1 (12) 
E V dS + 
én, ”(R(p;q) 
a l 
* fV a RoD 


Here, from the discretization of the body and wake, and 
assuming the potential @ and the value of 0@/On are constant 
within each panel. Then, equation (12) can be written by: 


2 {K M 
+2} 2, 2 (Ab, Vwi [+ (13) 
Np=1 | k=1_j=1 
k 
2 K Nig 
C 
+X /EE(S]ys5 
Np=1 k=1 j=1 on j 
— gees Niot 
Where V. A oa) ; Wi, and V, S* are the velocity influence 


coefficients. Those velocity coe Hicients can be evaluated 


analytically by assuming that the surface elements are 
approximated by a number of quadrilateral hyperboloidal 
panels. 

Calculations of the velocity influence were more sensitive 
than the potential coefficient, and also the required storage was 
three times more than the storage of the potential coefficient. 
There was one big advantage that the velocities can directly 
be obtained for any field points 

The induced velocity diagram of the CRP is shown in 
Fig. 2, where u, and u, denotes the axial and circumferential 
induced velocities, respectively. Since the interaction between 
two propellers, the total induced velocities may be expressed 
as follows: 


u); = (ui T (u); 
(14) 
(u), 7 ui E (u); 
For two fore and aft propellers, it is given: 
(u), = U) 3 0), = Uo F U) ae 


(u), = u), ; 0), = Gt Gs 


Where (u); implies the axial or circumferential velocities 
ati-th propeller induced by the j-th propeller. The subscript 1 and 
2 denote the forward or aft propeller, respectively. We observe 
is very small and neglected. 


that (u,),, is zero and (u), 


Wala 
Aft propeller 
Fig. 2. Relative Velocities at Blade section of CRP 


The hydrodynamic pitch angle and resultant velocities to 
the fore and aft propellers are expressed as follows: 


Vag +U, 
tan B= Adv aF 
OFf + Utp 
(16) 
tan B; 4 = 
Orp — Ua + (Uy )o1 


Inflow velocity to the fore propeller may be obtained by: 
Vir = Vaav + Uar + Op X Tp — Uig 


5 5 (17) 
Vir =VOWVaavt Uap) + 2apny 


a Urp) 
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Inflow velocity to the aft propeller is expressed by: 


(18) 


Via=Vadvt Uaa + (Ug)21 + @ 4 X Tap Ua + (Uy) 24 


2 3 
LA 5y (Vav t Uaat (Uaa) + Caran ao Urat (Uda) 


NUMERICAL RESULTS AND DISCUSSION 


Application of the CRP system 
on the bulk carrier and VLCC 


We applied the method on two CRP for two different ship 
types where the model tests have been done in IHI and were 
available [9] and [10]. Main dimensions for bulk carrier and 
VLCC are given in Tab. 1. 


Tab. 1. Main dimensions of the Bulk Carrier and VLCC 


Bulk Carrier 


Parameter Cuno) 


Length (LBP)[m] 


Breadth [m] 


Depth [m] 


Dead Weight 


Speed at full load [Knot] 


For each vessel, conventional propeller and CRP have 
been used and the results are compared. The CRP5022 and 
conventional propeller (single propeller MP588) for the bulk 
carrier and CRP5029 and MP620 for VLCC have been selected. 
The principal of particulars of the both CRP propeller and single 
propeller are given in Tab. 2 and Tab. 3. 


Tab. 2. Principal particulars of CRP and single propellers for Bulk Carrier 
i Contra-Rotating 
Propeller.| Single | Propeller CRP5022 
Type | Propeller 


D [mm] for model tset 


Parameters 


D [m] for ship 


Boss Ratio 


Skew angle [deg.] 


Direction of Rotation | Left hand | Left hand 


Section 
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Tab. 3. Principal particulars of CRP and single propellers for VLCC 


Contra-Rotating 


Propeller. Single Propeller CRP5029 


Type | Propeller 
MP620 


Parameters Foreward 


D [mm] (model) 


D [m] (VLCC) 


Boss Ratio 


Skew angle [Deg.] 


Direction of Rotation | Left hand | Left hand 


Section 


Grid generation 


As known that the BEM is dealing with boundary of the 
body, so the grid generation is the necessary starting point for 


Fig. 3. Surface Mesh on CRP5022 


I, 
f 
ra 


` 
Ati i 


Fig. 4. Surface Mesh on CRP 5029 


numerical implementation. The coordinate of the discretized 
surface should be sufficiently accurate since any inaccuracy can 
lead the hydrodynamic pressure to become noisy. According 
to our experience, the hyperboloidal quadrilateral element 
seems to be better than other elements like triangular and super 
element. Higher order element (quadratic or cubic order) is 
very complicated way to apply although it gives more precise 
results. This is our future plan to be done. However, for the 
present calculation and prepared SPD code, the hyperboloidal 
quadrilateral elements are used to discretize the whole body 
(hub and two propellers). 


HYDRODYNAMICS 
CHARACTERISTICS 
AND DETERMINATION 
OF DESIGN POINT 


Numerical results of the open water characteristics 
are compared with the experimental data for the CRP and 
conventional propellers. It is shown that the numerical 
results of the present method are very accurate and in good 
agreement with the experimental data for the open-water 
characteristics. 

Using the equations (1) and (2),(K,=A-J’; K,.=B-J’), in the 
figures of 5 and 6, intersection points are obtained between the 
thrust coefficients for each conventional and CRP propellers. 
Tab. 4 is shown the efficiency of each propeller (conventional 
and CRP) obtained from the figures 5 and 6 at the intersection 
points of Kt. 


—o— 10Kg-single propeller(Comp) 
i Kt-CRP(Comp) 
—e— Eta-CRP(Comp) 


== Kt-single propeller(Comp) 
4 Eta-single propeller(Comp) 
—*— 10Kq-CRP(Comp) 


59 
onesies oo = Kt-single propeller(Exp ) -<<- eee propeler Erp ) 
= = = = Eta-single propeller(Exp.) Kt-CRP(Exp.) 
aimeo 10Kq-CRP(Exp.) ---- Eta-CRP(Exp.) 
an Kt=A*I2 =W= Kt=B*J^2 


Kt, 10Kq. Eta 


Fig. 5. Comparison of open water characteristics of conventional and CRP 
propellers for the Bulk-Carrier and the determination of Design point 
~— Kt-single propeller(Comp) —a— 10Kq-single propeller(Comp) 
=a Eta-single propeller(Comp) = Kt-CRP(Comp) 
—s— 10Kq-CRP(Comp) —e— Eta-CRP(Comp) 


om =~ Kt-single propeller(Exp.) ------ 10Kq-single propeller(Exp. } 


MP620 ~~" Eta-single propeller(Exp. ) Kt-CRP(Exp.) 
CRP5029 7777 10Kq-CRP(Exp.) stents Eta-CRP(Exp.) 
nia =A= Kt=A*I*2 woe Kt=B*J^2 


Kt, 10Kq, Eta 


0 ee y y = = š = e : 
0.1 0.15 0.2 0.25 03 035 04 0.45 0.5 0.55 0.6 0.65 0.7 0.75 08 


J 


Fig. 6. Comparison of open water characteristics of conventional and CRP 
propellers for the VLCC and the determination of Design point 


Tab. 4. Comparison of Experimental 
and Calculated Efficiency for Types of Ships 


Ship Type Bulk-Carrier 


Propeller 


CRP5029 CRP5022 


J... = 0.62 


Exp 


Je, = 0.61 


0.63 


CONCLUSIONS 


In this paper, we numerically calculated the open water 
characteristics of the conventional and CCRP and obtained 
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the optimum operational condition for the large vessel using 

BEM. According to the results, following conclusions are 

drawn: 

O The present method can be applied to any complicated 
propeller configuration and determine the open water 
characteristics. 


O The CCRP system may raise the propeller efficiency around 
2-3 percent at design condition for the present system. 


O Design point is determined based on highest efficiency 
which is matched the generated propeller thrust and ship 
required thrust. 
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ABSTRACT 


The computer system presented in this article is composed of several program blocks for the complete design 
of ship propellers. The design calculations are based on a combination of the modified lifting line theory 
and on the vortex lifting surface theory. The system enables solution of the following design problems: 
- calculation of the scale effect on the ship wake velocity field, including the influence of the propeller 


and rudder on this field at the propeller location 


- maximization of the propulsive efficiency 


- optimization of the propeller blade geometry on the basis of the compromise between the cavitation 


and blade strength requirements 


- optimization of the number of propeller blades and blade geometry on the basis of the acceptable level 
of induced pressure pulses and unsteady shaft bearing forces 
- calculation of the blade spindle torque for the controllable pitch propellers. 


The computer system is equipped with many numerical options for graphical visualization of the input 
data, including an easy possibility of their correction and control of the intermediate and final results of 
calculations. 


Key words: ship propellers, design methods, computational fluid dynamics 


INTRODUCTION 


In the recent years a computer system has been developed, 
which enables the complete design calculation of the different 
ship propulsor types, including their analysis in the wake 
velocity field behind the full scale ship hull. For this purpose 
the design and analysis blocks of the system are supplemented 
with the block for calculation of the wake scale effect and the 
influence of the propeller operation and presence of the rudder 
on the ship model wake. 

The main design process is based on the modified lifting 
line model and on the model of vortex lifting surface theory. 
The analysis of the designed propeller operation in the non- 
uniform velocity field is based on the extensively modified 
program UNCA [2, 3, 4]. 

The propeller design system enables solution of the 
following design problems: 
© calculation of the effective ship wake velocity field, 

including the scale effect at the propeller location 
= correction of this field for the presence of the rudder 

(measurements of the wake in model scale are usually 

performed without the rudder) 
= maximization of the propeller efficiency 


= optimization of the propeller blade geometry, resulting 
from the compromise between the cavitation and strength 
requirements 

= optimization of the number and geometry of the propeller 
blades on the basis of the acceptable level of the pressure 
pulses generated on the hull and of the unsteady bearing 
forces 

= calculation of the blade spindle torque for the controllable 
pitch propellers. 


The system is based on the conversational principle and 
it is equipped with an extensive set of graphical programs. It 
enables an easy control and correction of the input data and of 
the results of calculations of every block of the system. All data 
modifications are recorded and they may be used, if necessary, 
in any further design calculation. 

The results of calculations may be presented in the form of 
printouts, graphical diagrams or films, showing the variation of 
the pressure distribution or cavitation phenomena on a rotating 
propeller. 

The system includes three main program blocks: 
> the program for determination of the design velocity field 
> the program for propeller design, 
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> the program for analysis of propeller operation in the non- 
uniform inflow field. 


The block diagram of the system is shown in Fig. 1. 
The expected final effect of the propeller design process 
may be achieved only through proper interaction between 
these programs. The calculations may be controlled from the 
computer screen, without preparation of the special data sets 
for the respective blocks of the system. The special graphical 
procedures enable an easy control of the input data and of the 
intermediate results, including also correction of the propeller 
blade geometry directly from the computer screen at all stages 
of the design process. 


Fig. 1. The block diagram of the computer system 
for the complete design of ship propellers 
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DESCRIPTION OF THE MAIN BLOCKS 
OF THE SYSTEM 


The input data 


The input data include all quantities necessary for the 
initiation of the four alternative versions of the calculations: 
* propeller design calculation without correction of the 

velocity field and without the analysis of the propeller 

operation in the non-uniform velocity field, 

* propeller design calculation without correction of the 
velocity field, but including the analysis of the propeller 
operation in the non-uniform velocity field, 

* propeller design calculation including the correction of the 
velocity field for the scale effect and the analysis of the 
propeller operation in the non-uniform velocity field, but 
without the correction for the presence of the rudder, 

* propeller design calculation including corrections of the 
velocity field both for scale effect and rudder presence 
and including the analysis of the propeller operation in the 
non-uniform velocity field. 


The input data may be introduced in the form of the pre- 
prepared input data file or they may be introduced directly 
from the computer screen in the conversational mode. The 
program is equipped with graphical procedures for control of 
the correctness and for modification of the input data. 


The program block for calculation of the scale 
effect and of the rudder and propeller influence 
on the wake velocity field 


In order to perform the calculation including the scale effect 
on the wake velocity field the appropriate data describing the 
hull geometry (typically the set of theoretical hull frames and 
outlines of the bow and stern) must be available. Alternatively, 
the earlier prepared data file containing the hull geometry 
defined as the set of panels may be used as input. The selection 
between these two options is made from the screen. In the first 
option the program for transformation of the hull geometry into 
the set of panels is used. 

Calculation of the scale effect on the wake velocity field 
does not require a large number of panels describing the hull 
geometry. The practical experience shows that the sufficiently 
accurate results may be obtained using: 

+ 60 panels along the hull length 
+ 20 panels along the hull frame (at one side of the hull). 


Such a number of panels may be accurately defined if 
the set of 16 theoretical frames is available, with more dense 
distribution in the bow and stern regions. Fig. 2 shows an 
example of the hull represented by the set of panels, viewed at 
a certain rotation angle. Such a presentation enables efficient 
control of the distribution of the panels, which should have the 
form of quadrangles possibly close to squares. 


Fig. 2. The ship hull viewed at a rotation angle around the x axis 


The correct representation of the hull as the set of panels 
allows further calculation using the PANSHIP program [9, 10, 
11, 12], in order to determine the scale effect on the hull wake 
velocity field at propeller location and to calculate the propeller 
design velocity field. 


The hull represented by the set of panels may be 
supplemented with the rudder described by a separate set 
of panels. The rudder is usually not taken into account in 
calculations for model scale (if the model experiments have 
been performed without the rudder), but it is always included 
in the full scale calculations. 


The determination of the design velocity field 


Taking into account the scale effect on the velocity field at 
propeller location is a very important aspect of the propeller 
design process. The design procedure itself, which determines 
the propeller performance and blade geometry, requires only 
the radial distribution of the circumferential average of the 
axial component of the velocity field. However, the complete 
design process requires information about the entire velocity 
field (three components distribution over the propeller disc) 
in full scale, including the effect of propeller operation (the 
effective field) and the effect of the rudder presence. This 
may be achieved only using the computation procedure 
based on the “computer model basin”. Such a procedure has 
been developed on the basis of the computer model basin 
PANSHIP [10, 12]. 

Example of application of such a procedure is presented 
in Figs. 3 and 4. Fig. 3 shows the results of measurements of 
the axial component of the velocity on the tanker model, while 
Fig. 4 shows the corresponding design velocity field calculated 
for the full scale ship. The influence of rudders of different 
thickness on the velocity field at propeller location is shown 
in Figs. 5 and 6 [6, 12]. 
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Fig. 3. Results of measurements of axial velocity component 
behind the tanker model Lpp = 6.53 [m] 


The results presented in these figures show the meaningful 
scale effect and visible influence of the rudder presence on the 
velocity field at propeller location. The analysis of theses results 
leads to the conclusion that either model experiments should be 
conducted with the rudder or the appropriate corrections may be 
calculated using the computer model basin (e.g. PANSHIP). As 
the velocity field should be corrected both for the scale effect 
and for the rudder effect, these corrections may be integrated 
in one procedure, leading to the effective full scale velocity 
field. Only such a velocity field enables correct analysis of the 
propeller operation behind the ship hull. 
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Fig. 4. Results of calculation of the design velocity field 
Jor the full scale tanker Lpp = 200 [m] 
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Fig. 5. Influence of the rudders of different thickness 
on the axial velocity component at the propeller location 
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Fig. 6. Influence of the rudders of different thickness on the radial 
distribution of the circumferential average of the axial component 
of velocity at propeller location 


The design program block 


The algorithm of the design program has been described in 
great detail in [1, 8]. The program based on this algorithm has 
been functioning over many years and it has been extensively 
verified. 

The design calculations are performed in all cases, 
irrespective of the version selected in the Section 2.1 above. 
In case of version a) these are the only and final calculations, 
ending with the results which are presented in the form of tables 
(cf. Figs. 7 and 8) and pictures shown on the computer screen 
(cf. Figs. 9, 10 and 11). The pictures may be easily printed. The 
full results are included in the file Design.OUT. 

In the case of version b) the design calculations are 
conducted alternating with the analysis program UNCA, until 
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the appropriate criteria regarding e.g. cavitation phenomena or 
pressure pulses and unsteady bearing forces are fulfilled. 


Fig. 9. The rendered picture of the designed propeller blades 


Fig. 7. The main results of the design task 
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0.3000 D. 0.1807 3.03003-01 
0.4000 0.3672 0.2002 -0.0166 
0.5000 0.3955 0.2141 -0.0163 4. 70002-02 
0.6000 0.4209 0.2169 ) -0.0085 
0.7000 0.4410 0.2155 0.0050 -9.00003-03 
ee oe oe oe > “+. socoz-oz 
1.0000 0.0430 -0.1285 a 4 0.1500 80.0000 -L.z1008-01 
Fig. 8. The table of the results of design calculations (the first page only) -4. 77008-01 
copied from the file Design.OUT -2.33008-01 
: -2.99008-01 
The program block for analysis of the propeller -3. 45008-01 
operation in the non-uniform velocity field ~4.01008-01 
. . -4. 57008-01 
, The program UNCA for the analysis of propeller operation =g.49008-43 
in the non-uniform velocity field is the very important element -5.69008-01 
of the propeller design process. The main part of the algorithm «6. 25008-01 
of this program deals with the determination of the extent and ~€.01008-01 
intensity of different forms of cavitation present on the blades -7. 37008-01 
of the propeller operating in the non-uniform velocity field. b M ~7.99008-01 
The original theoretical model integrates the unsteady vortex è pe ~6.<9008-63 
lifting surface theory with the dynamically varying cavitation can qeeeei 


bubbles. The detailed description of this program is included Fig. 11. The distribution of pressure coefficient on the suction side 
in [2, 3, 4]. of the blades in the propeller design condition 
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The input data to the analysis program are introduced from 
the appropriate file (the velocity field — either given or computed 
in the system according to Section 2.3) and from the results of the 
design task performed in Section 2.4 (geometry of the propeller 
and some additional data). These input data may be controlled 
from the screen as the Summary of the Analysis Task Settings. 
An example of such a table is shown in Fig. 12. 


+*ee* SETTINGS OF THE ANALYSIS TASK ***** 


t. THE BASIC PARAMETERS 
Selected values identical to the design task? 


Yes 
Technical parameters | 


ng the stern sou = (design) 

Inclination of the propelier s shaft [deg] = 

Off-design blade pitch setting a [eet = “<4 00 

Water kinematic viscosity [...] = 

Water surface tension [...] = 0. on 

Mean cavitation nuclei radius [...] 

Reference level for Re-calculation of the induced pressure pulses into decibels 
Technical parameters II 

Ship velocity [knots] = (design) 


+ SHIP MODEL VELOCITY FIELD FROM THE MODEL EXPERIMENTS 


- HULL OF THE FULL SCALE SHIP 
Model scale coefficient [-] = 28 
Block coefficient (initial) [-] = 0.6 
Location of the propeller relative to the hull 


. PRESSURES INDUCED BY THE PROPELLER IN SELECTED POINTS 
Calculation of the pressure pulses? 
Yes 


Fig. 12. Summary of the Analysis Task Settings 


The newly designed propeller is analysed using the program 

UNCA from the following points of view: 

%* presence of the different forms of cavitation in the selected 
angular propeller blade positions in the non-uniform 
velocity field 

%* values of the induced pressure pulses either on the ship hull 
or in the surrounding space, 

%* values of the unsteady bearing forces. 


After the analysis of the results of these calculations the 
appropriate modifications to the designed propeller geometry 
may be introduced and the design calculations may be repeated. 
For example the following propeller geometry parameters may 
be changed: 

+ the values and the radial distribution of the blade 
skewback 

the values and the radial distribution of the blade profile 
chord lengths 

the values and the radial distribution of the blade profile 
maximum thickness 

the type of chord-wise profile thickness distribution 

the type of chord-wise profile mean line camber 
distribution 

the radial hydrodynamic loading (circulation) distribution 
the number of blades. 


++ ++ + + 


The analysis of propeller operation in the non-uniform field 
of flow may be performed for the design ship speed and propeller 
rate of rotation or for the off-design values of these parameters 
(without changing the blade geometry and the velocity field). 
This is very convenient when the propeller is designed for the 
compromise design condition, e.g. for a fishing vessel, a tug 
boat or a navy ship. For such ships more than one ship speed is 
important (e.g. towing speed, free running speed or maximum 
speed) and the designed propeller should ensure the optimum 
performance over the entire range of operating conditions. 

An example of the selected results of the calculations of 
the unsteady bearing forces and cavitation extent is presented 
in the following figures. 


UNCAS3 LIFTING SURFACE PROGRAM FOR UNSTEADY PROPELLER CAVITATION ANALYSIS 


TKTAALTTTTTAAAAAT AKA TTSTL ALLA TTARAAAAAKSSSTLLLSKSSATAKAAAKAKE ERAT LLAe eee 


PROJECT NO= 1 
CALCULATED FOR=Destnax 
PAGE NO= 1 


CONTROL PRINTOUT OF THE INPUT DATA 


TKAAALTAATAAAAAATAKAESTSL LLL esses 


Kl= 1 Mi= 9 NZ = 5 {-) CVC = 1.1900*10-6 [N2S-1)} 
K2= 3 N2=15 VS = 11.3178 (MS-1) CP = 2500.00 (NN-2) 
K3= 1 m= 9 PD = 7.0000 (M) SDW = 1025.00 (KGH-3) 
K5= 3 N4=36 HD = 1.4000 (M) WTS = 0.0720 (NN-1) 
K7= 1 MS= 6 RPS= 1.5167 (S-1) RAVE= 0.0500%10-3 [N) 
KƏ" 2 N6=70 WSD= 4.5000 (M) RLP =0.100000 (PA) 
SFI" 0.0 (DEG) DLF = 0.0 (DEG) 
RR PI cr SI REKI THI CHI RLI 
0.20000 1.18110 2.13990 0.00000 0.00000 0.42630 0.12520 0.05000 
0.30000 1.19220 2.35970-0.08540 0.00000 0.32900 0.15570 0.04000 
0.40000 1.21420 2.57040-0.11620 0.00000 0.26180 0.17560 0.03500 
0.50000 1.21490 2.76650-0.11410 0.00000 0.21140 D.18440 0.03000 
0.60000 1.22080 2.94630-0.05950 0.00000 0.16800 0.18300 0.02500 
0.70000 1.26610 3.08700 0.03500 0.00000 0.12880 0.17260 0.02000 
0.80000 1.32650 3.15260 0.20510 0.00000 0.09170 0.15480 0.01500 
0.90000 1.40850 3.00930 0.48930 0.00000 0.05600 0.12880 0.01000 
1.00000 1.49040 0.30100 1.05000 0.00000 0.02800 0.01160 0.00500 


UNCAS3 LIFTING SURFACE PROGRAM FOR UNSTEADY PROPELLER CAVITATION ANALYSIS 


TERAATTAATAAAAAAAAHEEEEERES AA TATA AAAARKESERERAAS SAAT ANAK EERE ERRATA 


Fig. 13. The table of results of the analysis of the propeller in the non- 
uniform velocity field (the first page only) copied from the file Unca. OUT 


200 450 50 559 

Propeller position angie b [deg] 

Fig. 14. The components of the unsteady bearing force 
as the function of propeller angular position 


15.0 


100 150 20 20 %0 350 #0 450 S00 350 809 050 

Propeller position angle [deg] 

Fig. 15. The components of the unsteady bearing moment 
as the function of the propeller angular position 


Fx FY Fz MX MY MZ 


Fig. 16. The harmonic amplitudes of the three components 
of the bearing force and three components of the bearing moment 
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Apart from the static pictures of the calculated cavitation 
phenomena on the propeller (Figs. 17 and 18), the program 
can generate a moving picture of the cavitating propeller. In 
this film the propeller performs full revolution, showing the 
dependence of the cavitation phenomena on the local details of 
the inflow velocity field. All numerical results of the program 
UNCA are included in the file Unca.OUT. 


Sheet 
cavitation 


Bubble 
cavitation 


Intermittent 
cavitation 


Tip vortex 
cavitation 


Fig. 17. The example of visualization of the calculated 
cavitation phenomena on the propeller 


Sheet 
cavitation 


Bubble 
cavitation 


Intermittent 
cavitation 


Tip vortex 
cavitation 


Fig. 18. The frontal view of the calculated 
cavitation phenomena on the propeller 


FINAL REMARKS 


O The above presented system for the complete design of 
ship propellers enables an integral treatment of the hull- 
propeller-rudder system in the process of propeller design. 
The extensive graphical procedures for input data control 
and modification and for presentation of the results in 
the form of pictures and films are an important practical 
advantage of the system, facilitating its practical use by 
the designers. 


O The system has been developed with the support of the 
Research Grant No. 4T07C06630 of the Polish Ministry of 
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Science and Higher Education, together with other design 
systems for three special ship propulsor types, namely: 

+ ducted propellers 

* tandem co-rotating propellers 

* tandem contra-rotating propellers. 


O Description of these systems will be the subject of separate 
publications in the future. 
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Identification of shaft line alignment 
with insufficient data availability 


Lech Murawski, Ph. D. 
The Szewalski Institute of Fluid-Flow Machinery 
Polish Academy of Sciences 


ABSTRACT 


The paper presents a method of identification parameters of shaft line alignment in case 
of the lack of producers’ data. Proper shaft line alignment is often a problem for repair 
shipyards, for aged ships without sufficient documentation. Author proposed combined 
experimental-analytical method for identified some existing parameters and checking (and 
eventually correcting) power transmission system s foundation. Specialised software has 
%,| been developed for shaft line alignment calculations with influence coefficients. An example 
= analysis has been performed for cargo ships with medium-speed main engine and second 


one with slow-speed propulsion system. Multivariant computations supported by measurements of the 
ships’ shaft line have been carried out. 


Key words: shaft line alignment, bearings reaction, crankshaft external forces, shaft line stresses, aged ships 


INTRODUCTION 


Propulsion system’s foundation might be changed 
dangerously during ship hull repairing (especially welding 
works) even if there is no work with power transmission system. 
What is more, improvement of the shaft line - crankshaft 
alignment is difficult because very often aged ships haven’t 
available documentations. Repair shipyard does not know if 
foundation parameters of propulsion system are acceptable 
or not. Some of them (like shafts’ diameters, intermediate 
bearings’ reactions) are easy to measured but the others are 
practically inaccessible (shaft line deformation, shaft line 
— crankshaft interaction, stern tube and main engine bearings’ 
reactions, stresses). The analysis become much more difficult 
if our software take into account isolated power transmission 
system from ship hull (Fig. 1) without boundary conditions 
(e.g. without hull displacements and stiffness) [9]. Usually, 
similar parameters of shaft line alignment before and after 
repair is a shipyard’s target. 

There is no, well known, cheap methods for inspection 
and eventually improvement of propulsion system foundation 
after ship hull repairing [13]. Abnormalities of the propulsion 
system’s working parameters are often detected after ship’s 
repairing, during sea trial or even after certain period of 
ship’s voyage. Then, the failure repairing is very costly. 
Sometimes, tube of stern bearing or engine’s main bearing 
must be replaced and the proper shaft line alignment must 
be performed. 


Fig. 1. Typical model of the power transmission system 
isolated from ship hull 


POLISH MARITIME RESEARCH, No 1/2009 


35 


Author proposed a method of shaft line alignment 
identification based on easy to measured parameters and 
numerical analysis. Elaborated specialised software can 
calculate influence coefficients for each bearing. Shaft line 
alignment identification should be performed before and after 
repairing process. Analysis shows if correction is necessary (if 
shaft line alignment parameters have been changed significantly 
during repairing). If yes, analysis, based on influence 
coefficients, give us advice how easily improve propulsion 
system’s foundation. 


IDENTIFICATIONS’ METHOD 


Just before starting hull repairing some typical measurements 
[1, 6, 7] should be performed. First ofall, the jack-up test should 
be realized for intermediate bearings’ reaction identification. 
Also identification of shafts’ diameters is necessary. In the 
same time, the base of location of bearings body and shaft line 
position should be done. 

Position changes of the shaft line can be controlled 
(Fig. 2) by bending stresses measurement (e.g. by strain gauge 
technique) [2, 3]. Bending stresses should be recorded during 
full rotation of the shaft (with usage of turning gear) just before 
and after ship repair process. Changes of the stress level in the 
horizontal and vertical plane may be simulated by numerical 
analysis and then changes of the shaft line alignment may be 
identified. 


! HUS ae 
> Tnigemed titesbearine- 
~ Meastite ment placă 
r S n ° 
< E ofi shafi line 
bending Stresses 


Fig. 2. Bending stress measurement of the shaft line 


The relative position of bearing body might be based on 
measuring stress level [8] on foundation pads (Fig. 3) or on 
the other sensitive element of the bearing foundation (Fig. 4). 
Stress level should be recorded during jack-up test, before ship 
repair process. The simultaneous measurements can be treated 
as calibration process. After measurements the dependence 
between bearing reaction and stress level ofthe foundation pads 
(or foundation knee) is recognised. The same measurements 
(even without jack-up test) performed after ship repairing, give 
us changes of the bearing reaction which are induced during 
overhaul. 

Whole (before and after repairing process) measurements 
investigation gives us necessary data for shaft line alignment’s 
changes identification. The next step during expertises must be 
an assessment if identified changes are acceptable or not. The 
following limitations should be checked: loadings of the stern 
tube bearing, intermediate bearings and main bearings of the 
engine; stresses of the shaft line; interaction between shaft line 
and crankshaft. The shaft line alignment correction should be 
proposed if some limits are exceeded. Correction should be 
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practicable and cheap. Generally, it should be restricted to 
displacements of the intermediate bearings in horizontal and 
vertical plane. Author made a specialised computer program 
which can be useful for all this tasks (identification, assessment 
and correction) with limited (without design documentation) 
data availability. 


Fig. 4. Stress measurement of bearing s foundation knee 


SOFTWARE DESCRIPTION 


Assessment of the shaft line alignment for aged ships 
requires specialised software. The data for the software 
might be based on the measurements before and after ship’s 
repairing in case of insufficient design data availability. During 
ship hull repairing (especially during welding works) power 
transmission system’s axis might be unacceptable displaced 
even if there is no work with propulsion system. Assessment 
of the shaft line alignment before ship repairing process as 
well as determination of propulsion system’s sensitivity on 


stable hull deformation is very important from technical 
and economical point of view. Author performed specialised 
software for bearing reactions’ identification and shaft line 
alignment’s influence coefficients determination. 

Identification of shaft line alignment with knowledge 
only of bearings’ reaction and shafts’ geometry is reverse 
mathematical problem [4]. This problem is often ambiguous 
with multi solutions. For that reason equation no. 1 [5, 14] is 
solved during iteration process. 


([K] —[K,]) + {x} = tF} (1) 
where: 
K — global stiffness matrix, 
K,— geometric stiffness matrix, 
x — displacement vector, 
F — generalised force vector. 


Matrix equation no 1 is a system of algebraic equations with 
constant coefficients. Modified Gauss algorithm is used for 
solving the equations. The shaft line is modeled by 2-node (with 
6 degree of freedom for each node) linear beam elements with 
nonlinear boundary conditions. The beams can be cylindrical 
or conical and drilled or not. Bearings might be modeled by 
point wise or continuous support with taking into account 
bearing’s clearance and stiffness. Stiffness of the continuous 
bearing (e.g. stern tube bearing) is modeled by polynomial of 
10 order. In the program there are procedures for field solution’s 
searching. Finding out variants fulfilling measured values 
is a target of the algorithm. All variables are declared in the 
program as dynamic one so there are no limitations for number 
of elements and bearings. An example of multi bearings shaft 
line is presented on Fig. 5. 


ANALYSIS EXAMPLE FOR THE 
CONTAINER SHIP 


As a first step, the analysis has been performed for the 
container ship with typical propulsion system — constant 
pitch propeller driven directly by slow speed main engine. An 
example analysis has been made for container ship 4500 TEU. 
The hull stiffness characteristics on the base of FEM ship model 
analyses [10] like presented on Fig. 6. This kind of the model 
contains about 50 000 + 150 000 degrees of freedom. 

The methods of the characteristics determination of ship 
hull and main engine have been presented in the other author’s 
articles [11, 12]. The examples of ship hull with propulsion 
system dynamic characteristics are presented on Fig. 7-8. 

z|? 


1 9.000e+19 
0.000600 
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2.200004 


Fig. 6. FEM model of container ship 


Fig. 8. Transverse mode of main engine vibration 


Analysis is performed on the basis of standard shipyard’s 
shaft line alignment. This alignment is correct but their 
parameters might be better. The aim of the calculations is 
showing software possibilities. Standard shaft line reactions 
are presented on Fig. 9. 
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Legend of Dimensions: L- length: D- diameter: Da - aft end diameter: Df - fore end diameter. d - inner diameter 
Fig. 5. Shaft line model of two shafts universal supply ship 
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POZIE Z"9d 
COMIEST pug d 


INTERMEDIATE BEARING No. 1 


Reaction = 2.009e+02 

Shaft Deform. = 1.34001 mm 
Hull Deform. = 0.000e+00 mm 
Bearing Skew = 1.725000 mm 
Bearing Slacknass 


Reaotion Distribution 
Ratio = 0.00De+00% 


lacknes 5. 


INTERMEDIATE BEARING No. 2 


Reaction = 1 316e02 

Shaft Deform. = 2.774202 mm 
HA Detorn. = 0.000e+00 mm 
Bearing Skew = 4.273801 mm 
Bearing Siacknexs 

Ratio = 0.0000+00% 


for TUBE BEARING No. I t 


Ri= 3.47 leet? kN 

Ra= 6.465e402 KN 

Shaft Skew = -7.578¢-02 mm 
Hull Skew = 0.000400 men 
Resutant Skew = -7.578e02 mm 


Bearing Stackness 
NMaximum Ratio = 0 00000% 
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INTERMEDIATE BEARING No. 4 
INTERMEDIATE BEARING No.7 
Reaction = 2.484002 

Shart Deform. = 4.133-02 mm 
Hä Deform. = 0.000000 mm 
Bearing Skew = -6.631e-03 mm 
Bearing Stackness: 

Ratio = 0.000e+00% 


Reaction = 62230402 

Shaft Deto. = $.700e-02 mm 
Ha Derom. = 0.0008+00 mm 
Bearing Skew = 1.620002 mm 

Bearing Saokness 

Ratio = 0.000e+00% 


INTERMEDIATE BEARING No. 3 


Reaction = 2.742002 

Shatt Defo. = 1.628001 mm 
Hull Detorn. = 0 O00200 mm 
Bearing Skew = 34746-02 mm 
Bearing Stackness 

Ratio = 0.000e+00% 


INTERMEDIATE BEARING No. 5 
INTERMEDIATE BEARING No. & 

Reaction = 2,490e+02 

Shatt Detorn. = 4.103¢-02 mm 

Hull Deform. = 0000800 mm 

Bearing Skew = 2.307e-02 mm 

Bearing Stacknass: 

Ratio = 0.000e+00% 


Raaction = 6.112002 

Shaft Deform. = 8.506e-02 mm 
Wl Deform. = 0.000e+00 mm 
Baaring Skew = -3.621e-02 mm 
Bearing Stackness 

Ratio = 0.000e+00% 


INTERMEDIATE BEARING No. 6 


Reaction = 4.22702 INTERMEDIATE BEARING No. 8 


Shart Deform. = 7 0268-02 
Hull Detorn. = 0.000000 
Bearing Skew = 5.682002 mi 
Bearing Slackness 

Ratio = 0.000000 4 


Reaction = 2.0730+02 

Shaft Deform. = 3.460e-02 mm 
Hull Deform. = 0.000e+00 mm 
Bearing Skew = -1.068e-03 mm 
Bearing Slaokrass 
Ratio = 0.000e+00% 


Fig. 9. Container s shaft line reactions determined by the software for standard alignment 


The standard shaft line alignment might be improved 
by equalization of the aft stern tube bearing’s reaction and 
by reduction of the aft engine’s main bearing. Loading 
enlargement of the fore stern tube bearing (and aft stern tube 
bearing’s loading equalization) is a consequence of reaction 
reduction of the aft intermediate bearing (to 100 kN, ~24%). 
Resultant skew of the propeller shaft in the bearing tube is 
then decreased of about 8%. Loading reduction of the main 
bearings (of about 9%) is an effect of reaction increasing 
(to 300 kN) of fore intermediate bearing. The algorithm find out 
automatically new (for given intermediate bearings’ reactions) 
shaft line alignment. The modified parameters of the new shaft 
line alignment are presented on Fig. 10+12. 
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Fig. 10. Modified shaft line deformation 


Influence coefficients are also determined by the program. 
The coefficients are useful for designers because they show 
the relative changing of the shaft line parameters (shaft line 
deformation, reactions, bending moments, shear forces, 
stresses) if chosen bearing is moved vertically or longitudinally. 
An example of influence coefficients for vertical movement of 
aft intermediate bearing are presented on Fig. 13 + 15. Influence 
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coefficients of the bearings reaction under the influence 
of bearings’ longitudinal and vertical movements are also 


determined. 
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Fig. 11. Modified shaft line bending moments and shear forces 
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Fig. 12. Modified shaft line stresses distribution 
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Fig. 13. Influence coefficients of shaft line deformation under the influence 
of vertical movement of aft intermediate bearing 
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Fig. 14. Influence coefficients of bending moments and shear forces under 
the influence of vertical movement of aft intermediate bearing 


The analysed shaft line is relatively stiff — it is short with 
big diameter and the propulsion system is equipped in many 
lateral bearings. For that reason, displacement of one bearing 
have an influence on only one shaft segment limited by adjacent 


bearings. Longitudinal movement of the intermediate bearings 
has hardly perceptible influence on the analysed shaft line 
alignment parameters. 
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Fig. 15. Influence coefficients of shaft line stresses under the influence 
of vertical movement of aft intermediate bearing 


IDENTIFICATION OF THE MIDDLE-SPEED 
PROPULSION SYSTEM 


Identification has been performed for aged ship destined for 
heavy repair process. It is universal supply ship with double- 
shaft, middle-speed propulsion system (with gear box). The 
shaft line is more flexible then previous one: 282 is a diameter 
of intermediate shaft and 6400 is a diameter of propeller shaft. 
In the ship’s documentation, only design bearings’ reactions 
are available. Real reactions are known after measurements. 
Design and real shaft line alignment have been identified by the 
discussed software. Identified shaft line reactions are presented 
on Fig. 16 and 17. Significant difference between design and 
real reactions are observed. The most important differences 
(which needed improvements) are as follows: unloaded aft edge 
and overloaded fore part of stern tube bearing and insufficiently 
loading of stern intermediate bearing. Shaft line alignment 
should be improved. 

The author’s software determine given shaft line alignment 
as a first step, and searching the field of optimal solutions in 
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Fig. 16. Design shaft line alignment and reactions 
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Fig. 17. Real shaft line alignment and reactions 


the second step. In the third step the influence coefficients are 
determined by the program. The coefficients are useful for 
fast and precisely determination of proper bearings position 
in vertical and longitudinal plane. Determined influence 
coefficients are presented in Tab. 1. Denotations: ST — stern 
tube bearing, I1 — first (starting from the stern) intermediate 
bearing, etc. 


Tab. 1. Vertical influence coefficients 


Bearing: 


I2 [kN/mm] | 13 [kN/mm] | 14 [kN/mm] 


Target of the shaft line alignment is determining proper 
and uniform bearings reaction, acceptable stress level and 
good interaction between shaft line and crankshaft or gearbox 
[15]. Proper shaft line alignment has been determined (based 
on coefficients presented in tab. 1) by following lift of the 
intermediate bearings: I1 —4.8 mm, I2 —6.1 mm, I3 — 4.5 mm, 
I4 — 2.0 mm. Improved bearings reaction distribution is 
presented on Fig. 18. Comparison of the shaft line alignment’s 
parameters, before and after improvements, is shown on 
Fig. 19 + 21. 

Shaft line alignment after correction is characterised by 
smoother axis deformation and uniform bearings reaction. 
Intermediate bearing no. 1 is heaviest loaded by the others 
because of desirable loading of stern tube bearing. After 
correction loading of the fore edge of stern tube bearing 
is reduced two times, while loading of aft edge become 
acceptable. Before correction the aft edge of stern tube 
bearing was unloaded and was threaten by hammering 
phenomenon. 
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Fig. 18. Shaft line alignment and reactions after correction 
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Fig. 19. Shaft line deformations before and after correction 
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Fig. 21. Shaft line stresses distribution before and after correction 
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SUMMARY 


Changes of shaft line alignment of aged ship with 
insufficient data availability have been difficult for proper 
realisation with standard methodology. Shipyards are looking 
out for well experienced workers because shaft line alignment 
has been performed by tests and errors method. Usually, it is 
very costly and the cost ofthe process is depended on engineer’s 
experience. Sometimes, repair shipyards make an assumption 
that bearing reaction deviation + 50% is acceptable! 

In the paper author proposed a quite simple and cheap 
method for identification of shaft line alignment parameters 


in case of the lack of producers’ data. It is proposed combined 
experimental-analytical method for identification and 
correcting some existing parameters of power transmission 
system’s foundation. Specialised software has been developed 
and verified for shaft line alignment 
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Ring thruster — a preliminary optimisation study 


Leszek Matuszewski, Ph.D. 
Gdansk University of Technology 


ABSTRACT 


The ring thruster is a new type of propeller, for which there is no experimental data to 
verify analytical design calculations. A significant feature of the ring thruster is the absence 
of a shaft. Propeller blades are mounted to the ring rotating inside the housing, which 
has the shape of a nozzle. For this reason the ring thruster is closest, with respect to both 
the construction and principle of operation, to the Kort nozzle propeller. The absence of 
a shaft and no gap between the blades and the nozzle make it impossible to fully relay on 
results obtained from Kort nozzle propeller examination. What is more, the already existing 
SOPU codes developed for designing Kort nozzle propellers cannot be directly used for designing ring 
thrusters either. That is why for this purpose a new code determining hydrodynamic characteristics based 
on the theory of the vortex lifting surface will be used. When using the above method, some differences 
between calculated and experimentally recorded results are expected to be observed. To a significant 
extent, the level of the torque taken by the thruster will be affected by drag of the rotating ring to which 
the blades are fixed. Examining a propeller equipped with a rotating ring has revealed that the expected 
torque increment may reach as much as a few per cent, at the comparable level of axial force (thrust). 
At the present stage of ring propeller investigations there is no data available on how to shape the ring 
propeller blades. Possible comparison calculations, done using the existing computer code, will allow, the 
most, the shape of the blades to be determined for preliminary tests in the cavitation tunnel and on a self- 
propelled model. And only the results obtained in these tests will provide opportunities for verification of 
preliminary design calculations. It should be stressed, however, that developing design procedures for this 
type of propellers will require additional optimisation calculations, with further experimental verification. 
And this should be the subject of separate investigations. 


Key words: hydrodynamics, marine engine, thruster 


GOAL AND SCOPE OF EXAMINATION RESEARCH RIG 


The goal ofthe reported activities was to perform preliminary The research rig, which is the cavitation tunnel K11-MH 
hydrodynamic tests of the prototype ring thruster delivered by IMP PAN, was prepared for examining the new generation ring 
the employer. The main dimensions of the examined thruster thruster (a general view of the rig is shown in Photo 1). 
were the following : i 


* outer ring diameter Dzd = 0.208m 
* inner diameter of ring motor Dw = 0.130m 
* width of housing (nozzle) B = 0.1365 m 
* number of blades Z =6 


The scope of the examination includes measurement of the 
axial force induced by the ring thruster at different propeller 
revolutions and different velocities of the water approaching 
the thruster. 

A tensometric dynamometer JK-21-2-500N-2003 IMP PAN 
allows total thrust force of the Tpx unit to be measured. This 
force comprises the thrust generated by the propeller and the 
nozzle, and the drag generated by the elements connecting the 
thruster with the dynamometer (connectors and a cable). 

The torque is determined from the power calculated based 
on measurements of parameters of the electric current supplying 
the rotor, and therefore it can be burdened with some error. Photo 1. Research rig for ring thruster examination 
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A prototype strain-gauge dynamometer JK-21-500N-2003 
IMP PAN, mounted with relevant accessories and fixing system, 
is shown in Photo 3. 

The system for measuring, storing and processing of 
changing forces Tpx and Rox makes use of a microprocessor 
measuring amplifier AWO 100, linked with the computer via 
RS 232C interface. 

The average velocity in the measuring section of the 
cavitation tunnel was determined from the pressure difference 
measured in the tunnel confusor using a liquid-column gauge. 

The propulsion system of the ring thruster consists of the 
nozzle-shaped housing and the internal ring with six blades 
fixed to it, (Photo 2). 

Revolutions of the propulsion motor were controlled using 
a programmable inverter. 


Photo 3. Ring thruster in operation 


EXAMINATION PROCEDURES 


The substitute ring thruster (Dw = 0.13 m) delivered by 
the employer was mounted in a strain-gauge dynamometer 
JK-21-2-500N-2003 IMP PAN, especially manufactured for 
this purpose, and then installed on the properly instrumented 
research rig No. 2 (Cavitation Tunnel, type K11-MH) in the 
Szewalski Institute of Fluid-Flow Machinery PAS, Gdansk. 
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Ring thruster examination included the following tests: 

+ drag generated by a dummy ring thruster at different water 
velocities 

+ mechanical drag generated by the ring thruster in the air 
for different propeller revolutions 

+ ring thruster examination for different water velocities and 
different rotational speeds of the propeller. 


During the examination the following quantities were 

recorded: 

¢ Propeller revolutions n [rev/min] 

¢ Speed of water in tunnel measuring section V [m/s] 

¢ Total thrust and drag generated by ring thruster Tpx [N] 

+ Drag generated by dummy ring thruster, console, and cable 
Ro [N] 

+ Total torque taken by thruster Q [Nm] 

¢ Internal ring thruster drag torque in air Qo [Nm] 

+ Other recorded quantities included electric current intensity 
and voltage generated by the inverter supply system in the 
ring thruster motor. 


The examination was performed for steady water velocities 
and changing ring propeller revolutions. 

The examination was divided into test series, and the results 
of measurements were stored in the computer, in files named 
for instance ,,1200-01.dat” which means the first test series at 
revolutions 1200 rev/ min. 

Instantaneous torque values were stored in computer’s 
memory and then averaged, while the thrust generated by the 
thruster were averaged in the microprocessor based measuring 
instrument, by introducing a filter. The averaged results were 
then presented on a digital display. 


Photo 4. Research rig during ring thruster Siamination 
RESULTS OF EXAMINATION 


The tables below collect selected results of measurements, 
recorded on the research rig in the cavitation tunnel with the 
ring thruster. These results are grouped in properly numbered 
series. 

Quantities mentioned in the tables : 

Water velocity at the entrance to the measuring section is 
given in mm water column Ah. Real velocity is determined 
from the relation: 
=> V=0.1412 VAh [m/s] 
= total thrust (measured) Tpx generated by the entire set is 

given in [kG] 
= torque is given in [kGm] 
= revolutions are given in [ rev/min]. 


Tab. 1. Results of ring thruster examination 


F No. 
File L of 
[rev/min] 


series 


Tpx | Ah Q 
[kG] | [mm] | [kGm] 


1200-01 .dat. 1200 11.4 15 1 


1 
1 


800 

700 

1250 
1100-02.dat. 
1000-02.dat. 
800-02.dat. 
700-02.dat. 
1200-02.dat. 
1000-02.dat. 
1200-02.dat. 


1000 7.75 10 
1200 


1300-03 .dat. 
1200-03 .dat. 
1000-03 .dat. 
800-03.dat. 
700-03.dat. 


1300-04.dat. 
1200-04.dat. 
1000-04.dat. 


1400-05 .dat. 
1300-05.dat. 
1200-05 .dat. 
1000-05 .dat. 
800-05.dat. 


800-06.dat. 
1000-06.dat. 
1200-06.dat. 
1300-06.dat. 
1400-06.dat. 


1400-07.dat. 
1300-07.dat. 
1200-07.dat. 
1000-07.dat. 
800-07.dat. 
1500-07.dat. 


1450-08 .dat. 
1400-08 .dat. 
1300-08 .dat. 
1200-08 .dat. 
1000-08 .dat. 
800-08.dat. 
1477-08 .dat. 


Tab. 2. Torque generated by ring thruster working in air Q, 


n [rev/min] 


MOM 1400 
MOM 1200 
MOM 800 


The results of measurements, obtained after doing relevant 
recalculations to dimensionless coefficients of thrust, K, and 
torque, K are given in Figs 1 and 2. 

Fig. 1 shows the results of measurements of thrust (axial 
force) induced by the entire ring thruster, in the arrangement: 
dimensionless thrust coefficient K, vs. advance coefficient J. 

Fig. 2 presents the results of measurements of the torque 
taken by the ring thruster, in the similar arrangement: 
dimensionless torque coefficient K, vs. advance coefficient J. 
The dimensionless coefficients are defined as: 


K,= T/(pn’D*) 
K= Q/(pn?D’) 
J = V/D) 


Figs 1 and 2 show two curves, of which the broken curve 
refers to the thruster working in water. The continuous curve 
was obtained by deducting torque Q, , representing internal 
drag, from the total torque Q. 

1.6 


0 0.1 0.2 0.3 0.4 : 0.5 0.6 0.7 0.8 0.9 


Fig. 1. Results of calculations of thrust coefficient K, 
induced by the entire thruster 


0 0.1 0.2 0.3 0.4 : 0.5 0.6 0.7 0.8 0.9 


Fig. 2. Results of calculations of torque coefficient K ü 


The scope of measurements included thruster operation 
above J = 0.2 (measuring points are shown in the figures). For 
J lower than 0.2, including J = 0, the curves were extrapolated 
using the analogy to other Kort nozzle examinations. Using 
the same analogy the thrust induced by the entire thruster was 
hypothetically divided into part induced by propeller blades 
and that induced on the nozzle (ring). The proposed division 
is shown in Fig. 3. 
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0.5 0.6 0.7 0.8 


Fig. 3. Approximate division of total thrust K,, into part K, induced 
by thruster blades and part Ky induced by the surrounding ring (nozzle) 


The presented results should be treated as concerning 
preliminary investigations of this type of propellers in Poland. 
Very interesting results were obtained for total thrust and torque 
induced by the entire thruster. Diagrams in Figs 1 and 2 deliver 
a lot of new information on the performance of these types of 
propellers, but at the same time new problems can be addressed 
which need solving. 

What needs clarifying first is why the measured torque, 
induced on the blades and the ring, differs so dramatically (more 
than twice) from the calculated value of the torque induced 
by the blades alone. It would mean that the torque induced 
on the ring is very high and considerably affects thruster’s 
performance. If so, its reduction is a priority for improving 
thruster’s performance. 

It would also mean that the substitute seal and bearings of 
the examined thruster do not fulfil requirements concerning 
the minimisation of drag forces generated by the ring moving 
in the gap. 

In this type of investigations, good preparation and 
execution of measurement of a torque attributed to viscous 
losses on the ring, irrelevant of the torque induced by the blades, 
is the high priority. 

The analysis of thrust examination results (table 1) reveals 
that the ring surrounding the propeller has favourably affected 
the total thrust. Fig. 3 shows a diagram with hypothetical 
division of the total thrust K,, into part K,, induced by the 
blades and part Ki induced by the ring. 

The curve K,, was created using as a basis the point 
calculated for the blades alone and assuming the analogy to 
the Kort nozzle propeller examination. The shape of the curves 
is realistic, but it should be confirmed (or verified) by relevant 
examination of a propeller driven in a way classical for the 
Kort nozzle propeller system. 


CONCLUSIONS 
O To sum up, from the point of view of hydromechanics, 


independently of definite design solutions (including the 
motor, bearings, lubrication, seals) worked out for the 
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presented propulsion system, more comprehensive model 
investigations should be carried out to determine relevant 
empirical corrections which would allow these propellers 
to be designed in a way similar to that followed when 
designing Kort nozzle propellers. 

O High efficiency and dynamics of electric drive systems used 
on watercraft resulted in their increased proportion in total 
number of drive applications. Thanks to the development 
of mechatronics, electrotechnics and hydrodynamics, the 
time has come when earlier solutions in this area can be put 
in practice. Among water propulsion systems, especially 
attractive properties are represented by ring-type propellers, 
but only equipped with electromagnetic bearings. 
Acharacteristic feature of the motor, being an extension 
of a classical synchronous motor with permanent magnets, 
is that the ring with propeller blades is a part of the rotor. 
Dimensions of the nozzle in which the winding is mounted 
do not exceed dimensions resulting from the optimum 
geometry of the Kort nozzle propeller for the assumed 
power. As recently as a few years ago the development 
in the field of magnetic materials, ferrofluid liquids, non- 
linear control techniques and hydrodynamics reached 
a level providing opportunities for effective introduction 
of ring-type propulsion systems on a large scale. A drive 
which is expected to be especially promising in the nearest 
future is the ring thruster, used as a manoeuvring drive, 
and the main drive on smaller watercraft. The experience 
gained in designing ring thrusters with magnetic bearings 
has provides opportunities for manufacturing an efficiently 
working prototype and offering it on the market of water 
propulsion systems. 

O Ring propulsion systems with classical bearings are slowly 
becoming more and more popular, but still their efficiency 
is lower than 20% due to relatively high drag. A quantitative 
break-through in the field of efficiency improvement can 
be only secured by the use of magnetic bearings. 

O An additional advantage of the ring thruster is its ability 
to switch to turbine operation, in which it can be used 
as electric power generator. In numerous situations this 
property makes it possible to recover energy. 
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On the possible increasing of efficiency 
of ship power plant with the system combined 
of marine diesel engine, gas turbine 
and steam turbine, at the main engine - steam 
turbine mode of cooperation 


Marek Dzida, Assoc. Prof. 
Gdansk University of Technology 


ABSTRACT 


This paper presents a concept of a ship combined high-power system consisted of main 
piston engine and associated with it: gas power turbine and steam turbine subsystems, which 
make use of energy contained in exhaust gas from main piston engine. The combined system 
consisted of a piston combustion engine and an associated with it steam turbine subsystem, 
was considered. An algorithm and results of calculations of the particular subsystems, i.e. 
of piston combustion engine and steam turbine, are presented. Assumptions and limitations 
taken for calculations, as well as comparison of values of some parameters of the system 


and results of experimental investigations available from the literature sources, are also given. The system's 

energy optimization was performed from the thermodynamic point of view only. Any technical — economical 

analyses were not carried out. Numerical calculations were performed for a Wärtsilä slow-speed diesel 
engine of 52 MW output power. 


Keywords: combined systems, gas turbine, piston combustion engine, ship power plants, steam turbine 


INTRODUCTION 


For many years in domestic economy a trend to increase 
efficiency of both land and ship power systems has been 
developed. It is connected with economic problems (increasing 
fuel prices) and ecological ones (the lower fuel consumption 
the lower emission of noxious substances to atmosphere). 
New solutions are searched for to increase efficiency of power 
system by joining high-pressure combustion engines with 
other thermal engines, such as gas and steam turbines, into 
combined systems. 

In recent years combined power systems such as gas - steam 
turbine systems have been applied to ship power plants. Such 
systems in land solutions are capable of achieving the efficiency 
of more than 60 %. Ina ship solution such system was used on 
the passenger liner Millenium. However it has been so far an 
isolated example of the power system of the highest efficiency. 
The system requires to be supplied with a more expensive fuel 
oil, i.e. diesel oil (Marine Diesel Oil). Majority of merchant 
ships is driven by slow-speed diesel engines which combust 
arelatively inexpensive heavy fuel oil (Heavy Fuel Oil). It 
seems that such tendency will be still continued for many years 
in the world merchant fleet. 

High-pressure combustion engine (self-ignition one) is still 
the basic one used as a ship main engine, as on one hand it 
combusts the most inexpensive fuel oil (heavy oil), and on the 


other hand it has the highest efficiency as compared with all 
other heat engines. In exhaust gas from diesel engine is contained 
a huge amount of energy which can be utilized in another device 
(engine), that results in increasing the system’s efficiency and 
lowering emission of noxious substances to the atmosphere. 

One of such solutions represents a combined cycle of piston 
combustion engine and gas and steam turbines which utilize 
waste heat contained in exhaust gas from diesel engine. The 
piston combustion engine is the leading device in such system. 
It seems that now when fast containerships of 8 + 12 thousand 
TU capacity enter to service their propulsion systems require 
engines of a very large output power, of over 50 + 80 MW. On 
the other hand marine fuel oil is more and more expensive, and 
restrictive ecological limits for NO, and CO, emission to the 
atmosphere shall force to apply new power system’s solutions 
capable of increasing propulsion efficiency and lowering 
noxious gas emission to the atmosphere. 

For ship main propulsion large slow-speed marine diesel 
engines working on heavy fuel oil will be still used as their 
efficiency reaches today 48 + 51 %. At such large unit output 
values a huge amount of heat contained in exhaust gas and 
possible to be further utilized, is emitted from the engines. 

This paper presents a concept of a combined ship power 
plant fitted with a main slow-speed, piston combustion engine, 
as a leading one, and a steam turbine subsystem which utilizes 
waste heat from exhaust gas. 


POLISH MARITIME RESEARCH, No 1/2009 47 


A CONCEPT OF COMBINED 
SHIP POWER SYSTEM 


The assumed combined power system for large containership, 
having main slow-speed diesel engine and gas and steam 
turbines which utilize waste heat from the engine’s exhaust 
gas, is shown in Fig. 1. 

Exhaust gas discharged from particular cylinders of the 
main engine is gathered in the exhaust gas collector to supply 
the constant-pressure turbocharger. For the reason of high 
efficiency of turbochargers [1, 6], the power required for 
compression of supercharging air is achieved only from a part 
of exhaust gas flux. The remaining gas flux from the exhaust 
gas collector can be decompressed in an additional gas turbine, 
i.e. power turbine. In further considerations it was assumed that 
the power gas turbine drives additionally — through reduction 
gear — ship’s screw propeller. 

Exhaust gas from the turbocharger and power gas turbine 
flows to the waste heat boiler installed before the silencer in the 
main engine’s exhaust pipeline. The waste heat boiler produces 
steam used both for driving the steam turbine which delivers its 
power to the screw propeller, and for ship’s overall purposes. 


Fig. 1. A concept of combined ship power system 


Energy assessment of the combined 
ship power system 


The proposed concept of combined ship power system 
should be subjected to assessment of energy of its cycle. On 
the basis of the assumed schematic diagram shown Fig. 1, 
formulas for determining the power system’s efficiency, were 
elaborated. 

The power of the combined power system is determined by 
summing up power values of the system’s particular components 
(main engine, power gas turbine and steam turbine): 


N = ANT Na (1) 


combi 


hence the combined system’s efficiency is equal to: 
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and the specific fuel oil consumption is expressed by: 
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Np 

Np Pap - main engine’s efficiency and its specific fuel oil 

consumption, respectively. 
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From the relations (2) and (3) yields that every additional 
power in the power system makes the system’s efficiency 
higher and in consequence the specific fuel oil consumption 
lower. The higher power values achieved from utilization of 
waste heat contained in exhaust gas from main engine, the 
lower specific fuel oil consumption, therefore as large as 
possible output power values should be obtained from both 
the gas power turbine and steam turbine. Worth mentioning 
that ouput of gas power turbine depends to a large extent on 
turbocharger’s efficiency - hence a greater amount of exhaust 
gas can be decompressed in the gas power turbine (increasing 
its power) - as well as on efficiency of the gas power turbine 
itself. In steam turbine its power depends to a large extent 
on selection of live steam parameters as well as condenser 
parameters. 


Variants of combined 
ship power systems 


Depending on an assumed concept of combined system, 
systems consisted either of main engine and steam turbine, or 
main engine, gas turbine and steam turbine, can be considered. 
Contemporary design solutions of piston engine turbochargers 
require not very large amount of exhaust gas (due to their 
high efficiency). In such cases to apply the combined systems 
consisted of diesel engine, gas power turbine and steam 
turbine, is possible. This paper describes a combined power 
plant consisted of a piston combustion engine (serving as main 
engine) and a steam turbine subsystem. 


SELECTION OF MAIN ENGINE AND 
PRELIMINARY DESIGN ASSUMPTIONS 
FOR THE COMBINED SHIP POWER 
SYSTEM 


A 9RTA-96C Wärtsilä Sulzer slow-speed, two-stroke engine 
was assumed for calculations of the combined power system. 
The printout of its data calculated by means of the winGTD 
software [7] achieved from the engine’s producer, is presented 
in Tab. 1. 

For the so selected main engine the calculations of the 
combined system were carried out for tropical operation 
conditions on two main engine load levels (according to Tab. 1), 
i.e. 100% and 90 % CMCR (Contract Maximum Continuous 
Rating). 


Tab. 1. Parameters of 9RTA 96C engine [7] 


Np 
[%CMCR] 


be Mir 


[kg/h] 


[kg/kWh] 


379819 
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Fig. 2. The ship combined power system consisted of main slow-speed diesel engine and steam turbine (a concept) 


THE COMBINED SHIP POWER SYSTEM 
CONSISTING OF MAIN ENGINE AND 
STEAM TURBINE SUBSYSTEM 


The combined ship power system assumed for calculations 
in this section, consisted of the slow-speed piston engine and 
a steam turbine subsystem utilizing waste heat contained in 
exhaust gas from the engine, is presented in Fig. 2. Exhaust 
gas from the main engine turbocharger’s turbine flows into the 
waste heat boiler. The boiler produces the steam for driving the 
steam turbine which transmits its power to propeller shaft, as 
well the steam necessary for ship overall purposes. For the so 
assumed combined system, steam turbine subsystem (cycle) 
parameters were determined in two calculation points. 


CALCULATIONS 
OF STEAM TURBINE CYCLE 


A two-pressure waste heat boiler was assumed according to 
[3],[ 7]. Fig. 3 presents the schematic diagram of steam turbine 
subsystem with waste heat boiler. Exhaust gas flows from the 
turbocharger through the waste heat boiler transferring the heat 
to the steam cycle. In the high-pressure cycle (p,,) is produced 
the overheated live steam of the parameters (p,, t,) and mass 
flux m, which supplies the steam turbine. The low-pressure 
cycle (p) produces the saturated steam whose one part of m, in 
amount is extracted for ship overall purposes, and the remaining 
part is delivered to the steam turbine increasing its power. The 
applied steam turbine of condensating type additionally drives 
— through a reduction gear — the ship screw propeller. The 
assumed deaerating heater of boiling type supplies the water 
of the temperature t,,,, to the waste heat boiler. The condensate 
m, from the condenser and the condensate m, from the cooling 
box of the subsystem for overall ship purposes flows into the 
deaerating heater. An additional heat necessary to increase 
water temperature up to its boiling point in the deaerating 
heater is taken from the steam turbine’s extraction point; the 
heat is contained in the extracted steam flux m,,. The exhaust 
gas leaves the waste heat boiler at the temperature t, 

In each case it was assumed that the overheated live steam 
temperature at outlet from the boiler is lower than that of 
flowing-in exhaust gas by At, i.e.: 


t =t At (4) 


o inlet_B T 


The calculations were carried out in function of the pressure 
p, (equal to the p,, cycle pressure enlarged by hydraulic 
losses). 


Pa~ 1-05 p, (5) 


Fig. 3. The schematic diagram of the steam turbine subsystem with waste 
heat boiler. Note: for explanation of symbols — see Nomenclature 


Limitations of parameters and input values 
assumed for steam turbine cycle calculations 


Limitations of values of steam cycle parameters result 
from strength, technical and durability conditions of particular 
elements of the subsystem as well as from design and 
economical factors. The calculations of the steam subsystem 
of the ship combined power system were performed under the 
assumptions and limitations given below. 
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The temperature difference between exhaust gas temperature 
and that of live steam, At, for waste heat boilers of marine 
application, was assumed At = 15°C in compliance with 
the literature sources [5, 6]. The value of the ,,pitch point” 
recommended by the firm MAN B&W [8] for ship boilers is 
equal to dt = dt, = dt, ~ 10°C. The limit steam dryness ratio 
behind the steam turbine was assumed X, = 0.88. For ship 
outboard — water — cooled condensers the firm MAN B&W [8] 
recommends to assume the pressure inside the condenser equal 
to p, = 0.065 bar. Temperature of water supplying the boiler 
should not be lower than 120°C, at the sulphur content greater 
than 2% [8]. The reason is that outer surface of heater pipes 
from the exhaust gas side has its temperature higher by 8 + 
15°C than that of the supply water, and that materials of a higher 
resistance against acid corrosion are applied; hence for the 
calculations it was assumed that the supply water temperature 
can not be lower than t,,, > 120°C. It was also assumed that 
the exhaust gas temperature at outlet from the boiler must 
be higher by 15°C than that of supply water, i.e.: tp > tpw + 
+ 15°C. Every ship combusting the heavy oil for propulsion uses 
a large amount of heat to prepare the fuel. For the considered 
containership’s power plant the mass steam flow rate for power 
plant use in compliance with the MAN B&W recomendations 
[8] was assumed. The flow rate was made greater by the 
steam flow rate necessary for ship living purposes, which was 
estimated to be m, = 2000 kg/h of saturated steam. According 
to [8] the steam pressure for living purposes should be kept 
within the range p, = 7 + 9 bar. In the assumed scheme of steam 
turbine cycle, the pressure of the steam for living purposes is 
equal to that of the boiler’s low-pressure cycle: p, = p,.. The 
temperature in the cooling box was assumed t = 50°C, Fig. 3, 
in compliance with the recommendation [8]. 


Steam cycle optimization 


On the basis of the presented algorithm of steam turbine 
cycle calculations a computer software for determining the 
cycle parameters with the assumed quantities and limitations, 
was elaborated. Numerical calculations were carried out in 
variants, i.e. for: 


* the assumed pressure range of high-pressure cycle (of live 
steam): P, € P, min? P; a 

* the assumed pressure range of the boiler’s low-pressure 
cycle: Pi € (P, min? Pi au 

* the assumed pressure in the deaerating heater: 


p g € (p g min’ p baa 


Optimization of the steam cycle should be so performed as 
to make use of the heat contained in exhaust gas. Hence it is 
equivalent to determination of such steam cycle parameters at 
which the steam turbine reaches its maximum output power. 
The area of searching for optimum steam cycle parameters, 
shown in Fig. 4, should be narrowed down to that in which the 
above assumed limitations imposed on the steam subsystem 
are satisfied. 

The application of the steam subsystem with the waste 
heat boiler makes output power of the entire power system 
increasing both at 100% and 90 % CMCR load of the main 
engine. 

In the Tab. 2 the optimum steam cycle parameters for the 
two load levels of the high-pressure combustion engine, are 
presented. The maximum output power of the steam turbine 
obtained from optimization calculations, equal to N rma = 
= 3414 kW for the diesel engine load of 90%, is marked A 
in Fig. 4. No limitations were taken into account in this case. 
However in the real system limitations imposed on the steam 
subsystem should be taken into account. For the assumed 
limitations, i.e. the narrowed down area of permissible 
parameters of steam subsystem, the point of the maximum 
steam cycle output power N grma = 3182 kW, is marked o. The 
addition of the steam turbine to the high-pressure combustion 
engine system resulted in the increase of the power system’s 
output by AN N, = +7.43% and 6.87% for 100% and 90% 
load, respectively, at the simultaneous decrease of fuel oil 
consumption in the combined system from 174g/kWh to 
162 g/kWh, and from 169.8 g/kWh to 158.9 g/kWh, 
respectively, see Tab. 3. The specific fuel oil flux is lowered 
by 6.91% for 100 % load and by 6.43 % for 90 % load of 
the marine diesel engine as compared to those of the marine 
diesel engine alone. In the combined system its efficiency 


The optimum point without limitations 


p, [bar] 


Fig. 4. The permissible area of parameters of the calculated steam cycle 
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Tab. 2. The optimum steam cycle parameters in the combined system of high-pressure combustion engine and steam turbine 


P, tan 
[bar] | [PC] 


t 


FW 


[°C] 


without limitations 


ith limitations 


increased from 49.66% to 53.06% (An omb Mp = +7.43 %) 
for 90 % load and from 48.46% to 52.05% (Anom Mp = 
= +6.87%) for 100 % load. 

In the steam subsystem with limitations, for e.g. 90% load, 
the mass flow rate m, = 15.02 t/h of live overheated steam of 
the parameters: t = 289°C and p, = 19 bar, and m, = 7.15 t/h of 
live wet steam of the temperature t, = 165°C and the pressure 
p, = 7 bar, is produced from the heat contained in exhaust gas 
discharged from the piston engine, see Tab. 2. 


THERMODYNAMIC ANALYSIS OF THE 
COMBINED POWER SYSTEM 


The application of the combined power system to ship 
propulsion results in the increasing of propulsion system’s 
efficiency and — in consequence — in the lowering of specific 
fuel oil consumption and — additionally — in the increasing of 
output power without any additional fuel oil consumption. 
In Tab. 3 the obtained values of output power, efficiency and 
fuel oil consumption for the considered ship combined power 
system with limitations taken into account, are presented. 


Tab. 3. Output power, efficiency and specific fuel oil consumption 
of the combined power system consisted of 9RTA 96C two-stroke 
diesel engine and steam turbine 


ENGINE STEAM COMBINED 
RTA96C TURBINE | POWER SYSTEM 


AN /Np 


e combi ecombi eD 


[g/kWh]|  [%] 


ME AN oni Mp 
[%] [%] 


CONCLUSIONS 


O The proposed concept of the ship combined power plant 
consisted of slow-speed piston engine and additional steam 
turbine which uses the heat contained in exhaust gas from 
the piston engine, makes it possible to lower its specific 
fuel oil consumption by about 6.4 %. 

O Additionally, without any change of main engine power 
output the whole power plant output can be increased by 
about 7% due to utilization of the heat contained in exhaust 
gas. 

O 9RTA96C diesel engine consumes 205 t/day of fuel oil at 
96% CMCR power output, and in the case of application 
of the combined power plant the same power output is 
obtained at 90% CMCR load of main engine and the specific 
fuel oil consumption of 189 t/day. This way the shipowner 
obtains measurable profits in the lowering of the fuel oil 
consumption by 16 t/day, i.e. 7.8% savings in the daily fuel 
oil consumption. 

O Application of the additional steam turbine subsystem 
makes it possible to utilize waste heat from exhaust gas 
and in consequence to increase ship power plant output 
at a lower fuel oil consumption, and a lower emission of 
NO „ CO, and SO, compounds to the atmosphere without 
any additonal efforts. 


NOMENCLATURE 


- specific fuel oil consumption 

- specific heat of exhaust gas and air, respectively 

- specific enthalpy 

- mass flux of a medium 

- power 

- pressure 

- temperature 

- calorific value of fuel oil 

- efficiency 

- isentropic exponent of exhaust gas and air, respectively 


Indices 

bar - barometric conditions 

B - boiler 

combi - combined system 

D - marine diesel engine, supercharging 

f - fuel 

inlet - inlet passage 

k - parameters in a condenser 

o - live steam, calculation point 

air - air 

ss - ship living purposes 

T - stage of: compression in a compressor, decompression 
in a turbine 

C - compressor 

g - exhaust gas 
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T - turbine 

TC - turbocharger 

PT - power turbine 

ST - steam turbine 

exh - exhaust passage 

FW - water supplying a waste heat boiler 
BIBLIOGRAPHY 


1. 


2; 


52 


A new turbocharger generation. ABB Turbo Systems Ltd. 
Publication No. CH-Z 2035 98E 

Description of Marine GenSet- Type L28/32H. MAN B&W 
Diesel A/S. November 1993 


. Dzida M., Domachowski Z.: A combined gas-steam system in 


ship electric power system — prospects for large passenger liners 
(in Polish). In: [Materials of] 8th International Symposium on 
Ship Power Plants, Szczecin Maritime Academy, Szczecin 1997 


. ISO 8L17: Net specific energy of marine fuels. 1987 


POLISH MARITIME RESEARCH, No 1/2009 


. Kehlhofer R.: Combined-Cycle Gas & Steam Turbine Power 


Plants. The Fairmont Press, INC. USA 1991 


. Schrott K. H.: The New Generation of MAN B&W 


Turbochargers. MAN B&W Publication S.A. 236 5581E. 1995 


. Sulzer RTA 96C. Engine Selection and Project Manual. Wärtsilä, 


June 2001 


. The MC Engine. Exhaust Gas Date. Waste Heat Recovery 


System. Total Economy. Publication of MAN B&W. October 
1985. 


CONTACT WITH THE AUTHOR 


Marek Dzida, Assoc. Prof. 
Faculty of Ocean Engineering 
and Ship Technology 
Gdansk University of Technology 
Narutowicza 11/12 
80-233 Gdansk, POLAND 
e-mail: dzida@pg.gda.pl 


Photo: Cezary Spigarski 


POLISH MARITIME RESEARCH 1(59) 2009 Vol 16; pp. 53-60 
10.2478/v10012-008-0011-y 


Diagnosing turbine engine bearing structures 
with the aid of FAM-C and FDM-A methods 


Andrzej Gebura, Ph. D. 
Air Force Institute of Technology 


ABSTRACT 


The diagnostic method FDM-A [2], based on the measurement of the pulsation component 
frequency modulation, and the FAM-C method, making use of the alternating current 
frequency, were developed in the Air Force Institute of Technology. The essence of their 
accuracy consists in “natural” synchronisation of the sampling signal with the angular 
velocity of the observed kinematic links. These methods make it possible to determine the 
level of wear of a subsystem and its location during normal operation of the power unit. They 
provide opportunities for determining numerous parameters of bearings, bearing structures, 


and the rotor unit. Moreover, the FDM-A method makes it possible to detect resonances in particular elements 
and their parameters - among others, they provide opportunities for assessing the level of the mechanic 
quality factor of the examined structure, i.e. the time margin of structure operation until the structure 
reaches the state of resonance. On the other hand, the FAM-C method provides opportunities for defining 
numerous parameters of the phenomena connected with the rotor unit. The diagnostic system is connected 
to an arbitrary point of the electric network (electric junction, for instance) supplied with the direct voltage! 
or the alternating voltage? taken from the generators driven by the examined rotating system. 
The presented experimental material was recorded in operating investigations of the aircraft turbine engine 
unit, in which the role of “generators — observers” was played by: an aircraft DC generator and a three- 
phase AC rate generator. Generators of this type are also used in rocket cutters with turbine engines, this 
way the conclusions presented in the article can be easily transferred to those objects. 
During the measurements, performed on the decks of the examined airplanes and on engine test beds, 
the ranges and locations of assembly and wear defects were determined. Successive stages of wear of the 
rolling bearings were also observed. 


Key words: technical diagnostics, diagnosing propulsion transmission ducts, skew’, 
frequency modulation, characteristic set, turbine engine, rolling bearing 


INTRODUCTION 


Industrial development observed on the turn of the ninetieth 
century was a source of demand for rolling bearings, to be used 
in machines built in large numbers in those days. At the same 
time, the development of technology made it possible produce, 
in a mass scale, precisely machined bearing components. The 
bearing structures had relatively low rotational speeds, which 
rarely reached the level of 3000 rev/min. 

New challenges for rolling bearings were brought by the 
development of aviation — increased vibrations and the action 
of higher and rapidly changing temperatures made it necessary 
to secure higher reliability of the bearings via increasing radial 
clearances and more careful selection of construction materials. 
Even bigger challenge was connected with new constructions of 
jet turbine engines built after the second world war. Rotor units in 
these engines reached rotational speeds ranging from 10000 rev/ 
min to 60000 rev/min, much higher than ever before. 


are in use in the form of automatic tester). 
Also referred to as sweep. 


For a relatively long time, diagnosing the technical state 
of rotor unit bearing supports in aircraft engines has been 
done using vibroacoustic methods. The vibroacoustic signal, 
which carries the diagnostic information, is taken from the 
bearing structure elements and therefore is heavily noised and 
damped, which makes the diagnosis extremely difficult and 
time consuming. This is because of the presence of strong 
disturbances [5] generated by vibrating metal elements of 
the aircraft engine construction (which is, as a rule, relatively 
elastic). The second factor is the attenuation of the diagnostic 
signal on joints between the mechanical elements and in all 
spaces filled with air, fuel, oils and lubricants (all these are the 
media revealing poor propagation of acoustic waves). And the 
third, last factor is the concentration of numerous vibrating 
elements located close to each other, which makes detecting 
the location of the defected objects extremely difficult. 

In the article, another diagnostic method is proposed 
in which the diagnostic signal from the bearing supports is 


This method is called FDM-A (FM- frequency modulation, D - direct current, A - method advancement level) 
This method is called FAM-C (FM- frequency modulation, D - alternating current C - method advancement level- means that applications 
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transmitted via mechanical joints to an aircraft generator, in 
which it is converted to a synchronous, digital signal. This 
signal can be easily processes using present-day computers. 
Two variants of the proposed method were named FAM-C [2] 
and FDM-A by the authors. In those methods the analysed 
diagnostic signal is coded in the frequency and phase 
modulation, which is (as is well known from the theory of 
radio wave transmission) more resistant to disturbances than 
the amplitude modulation. The acoustic signal is subject to 
two-stage sampling: 
ve “synchronous sampling”, which consists in the conversion 
into a digital signal (instantaneous frequency time-history) 
by the “regular” deck generator 
ve “classical” sampling in time intervals between the above 
described “synchronous” sampling, done by a clock with 
properly selected frequency in the counter card. 


In Ref. [7-10, 13] the author has presented his achievements 
concerning the application of the FAM-C method making 
use of a single phase output voltage in the deck generator for 
diagnosing simple cases of wear and assembly defects of power 
transmission boxes. These achievements included: 
< monitoring of slips in one-way clutches, 
<+ estimating skew angles of power transmission shafts and 

their parallel displacements, 
~ detecting cavitation phenomena in the deck generator 

hydraulic rotational speed stabiliser system. 


Ref. [11] describes the use of a similar diagnostic method 
for estimating skew angles of power transmission shafts and 
their parallel displacements in the power transmission system in 
the ship engine room. The use of the three-phase measurement 
[13] made it possible to trace, without instrumentation 
change, mechanical processes revealing higher frequencies, 
for instance monitoring the operation of ship compression- 
ignition engines. 

Further resolution refinement was done using the signal of 
the DC generator output voltage. This method, named FDM-A, 
makes it possible to trace mechanical processes [12] revealing 
even higher mechanical frequencies: 

x monitoring the wear of toothed wheels, and estimating 
dimensions of tooth-to-tooth clearances 
X estimating the operation of rolling bearings. 


In author’s opinion, the here presented methods: FAM-C and 
FDM-A can become a valuable version of the non-interfering 
monitoring of the technical state of turbine engines. Some 
applications discussed in further sections indicate important 
advantages of these methods, such as: 
> no need for connecting any sensors to the object — this role 

is played here by the “regular” generator mounted on the 

object 

> the measurement can be done at an arbitrary network point, 
even extremely distant from the examined object 

> the primary diagnostic signal is directly “transmitted” from 
the examined bearing situated deep inside the structure to 
the generator serving as a sensor-processor via elements of 
the transmission construction. This way large disturbances, 
which can come from casing resonances, volumetric 
damping and other disturbances of high importance for 

vibroacoustic methods, can be omitted [5] 
> transient loss of the “connection” between the diagnosed 

element and the generator-receiver is not a problem — the 

clearances are presented as a separate estimate (among other 
actions, a separate characteristic set is created), which can 
also be an object of parametric analysis 
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> easy automation of the diagnostic process 
> high speed of diagnostic data recording. 


METHOD DESCRIPTION 


A comprehensive description of the diagnostic methods is 
given in Ref. [2,6 + 13]. It is advisable, however, to present its 
main elements, which are indispensable for understanding the 
context of the problems discussed in the article. Each assembly 
or wear defect, for instance skew of the splined connections in 
a power transmission box, is accompanied by the modulation 
of the output angular velocity. The modulation interval is 
a characteristic parameter for the type of defect and nominal 
angular velocity of the given kinematic pair. At the same time 
the level of frequency deviation (amplitude of modulations) 
is proportional to the scale of the defect. These modulations 
are transmitted, via the transmission system, to the rotor of the 
generator. Being a synchronous machine, the aircraft generator 
reflects changes of instantaneous angular velocity in the form 
of output voltage frequency modulation. By measuring time 
increments between successive zero level crosses and drawing 
their doubled inverses on the (t, f) rectangular coordinate 
plane we obtain a set reflecting, in a discrete form, the course 
of generator rotor velocity changes. Two parameters can be 
attributed to each deviation, which are: 
> deviation time t y; 
> deviation amplitude AF. 


The deviation time can be substituted by the process 
frequency Ís according to the formula: 


f = Mta) (1) 


Each deviation can be presented on the rectangular 
coordinate plane as a point (f,, AF). It was found in numerous 
research activities that these points tend to concentrate into 
sets. These sets bear the name of characteristic sets, as they 
characterise the wear of particular subsystems. It was found 
that they have different shapes, heights and locations with 
respect to the X-axis. 

It was also found that the change in the scale of the defect 
results in the change of the height {/AF,, / + /AF_,/}of the 
corresponding set, while the band 0-f, occupied by the set 
with respect to the X-axis is characteristic for the type of the 
subsystem. 

During numerous applications of the method the author 
has found that the presentations created in the above way as 
characteristic sets are favourable for a well known object, with 
well recognised relations between the change in the scale of 
the mechanical defect and the change of the characteristic 
set height. The use of these presentations has numerous 
advantages: 
= demodulation processes take a natural course directly 

in the deck generator, which simultaneously is a natural 

band-pass filter in the process of observation of mechanical 
phenomena 
= the electric signal recorded from the deck generator is 
synchronised in a natural way with the observed mechanical 
phenomena, as the rotor of the generator-observer is fixed 
to the diagnosed mechanical kinematic pairs 
= easy automation of the diagnostic process: 
¢ each type of defect of each subsystem has reserved an 
individual, earlier defined characteristic band on the 
X-axis 
¢ for each band, a number of diagnostic levels (usually 
four: A, B, C, D) are defined on the Y-axis — depending 
on the height of the given set the subsystem is classified 
as belonging to one of these diagnostic classes 


= easy observation of the time-histories of arbitrarily long 
duration in the same plane — of special importance when 
the signals appear stochastically 
= easy extraction of the pulsation component from a heavily 
modulated signal. Indeed, these extraction can also be done 
using partial Fourier analysis, but that analysis, unlike the 
FAM-C and FDM-A methods, has two disadvantages: 
¢ spectral lines for real power transmission systems are 
usually broadened, which makes the diagnostic pattern 
obscure and difficult for automation 
¢ that analysis neglects stochastic processes. 


Despite all above named advantages of the wear presentation 
on the (f> AF) plane, for unknown objects* it is advisable to 
use relations? f = f(t) for diagnostic purposes. Based on these 
relations, an experienced diagnostician can perform the process 
of state identification of the diagnosed object in a more efficient, 
although also more time consuming way. In that way stochastic 
processes can be traced, including, for instance, a temporary 
stop of the rolling element in its motion around the raceway 
due to the local appearance of the abrasion product. 

Once the characteristic sets are created based on the 
recorded time-histories, Fig. 3, the following parameters can 
be calculated: 

* height of the first harmonic characteristic set — proportional 
to the dimension of the total radial clearance and the rotor 
unit unbalance 

* width of the characteristic set — proportional to the rotor 
unit unbalance 

* distance between the first harmonic characteristic subsets 
— proportional to the shape error (ovalisation, roundness 
deviation) of the bearing seat 

* height of the main shaft first harmonic characteristic set 
— proportional to the rotating unit symmetry axis skew angle 
(between the generator shaft and the turbine shaft) 

* skews and axial clearances of other rotating engine 
elements, such as pumps, power transmission shafts, etc. 


DESCRIPTION OF THE DIAGNOSED 
OBJECT AND THE MEASURING SYSTEM 


Turbine engines have found their immovable place in 
aviation. An element which is most exposed to the action of 
dynamic forces in those engines is the bearing cage. According 
to the design tradition, its only role is to separate the rolling 
elements, as a result of which the bearing cage is not expected to 
carry heavy loads. The situation of general underloading of the 
bearing and rapid changes increasing its load, when the bearing 
cage is subject to the action of impact forces, is untypical for 
machine design traditions. Different kinematics of the rolling 
elements, being the result of different friction conditions, is the 
origin of impact forces in the rolling unit. The energy passed to 
the cage during its collision with the rolling element changes 
its momentum in both the circumferential and radial direction, 
and is partially lost for an additional slip of the rolling elements 
on their raceways. As a result of these phenomena, the cage 
starts to oscillate in an undamped way between the rolling 
elements [6]. It also rotates, in general, around another axis 
than the bearing axis. The stability of the cage motion highly 
depends on friction conditions observed in the rolling contacts 
and in the areas where the cage contacts the rolling elements. 
According to author’s experience, these forces are so large that 


4 
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they can provoke cracks of the separating elements in the cage, 
or even cracking of the bearing perimeter. Such phenomena are 
frequently preceded by resonance effects. 

Phenomena of this type, untypical for earlier applications 
of the rolling bearings, were the reason why producers and 
institutions involved in aviation industry often carry out, on 
a wide scale, research activities oriented on detecting and 
monitoring complicated processes of rolling bearing failures in 
aircraft engines. Nowadays, this problem is still unsolved and 
is an object of concern of designers, operators, and analysers 
of aviation accidents. 

In order to gain more detailed information on the relations 
between frequency modulation phenomena and bearing defects, 
the process of long-lasting wear of engine bearings was 
examined on the engine test bed, and was complemented by 
observations of engines directly on the decks of the airplanes. 
The object of investigations was a single-shaft turbine engine, 
Fig. 1. The main shaft is divided into the compressor shaft 
and the turbine shaft, which are linked together via the splined 
connection, Fig. 2. The shaft is supported on three rolling 
bearings. 


Fig. 1. Kinematic diagram of a one-watt turbine engine : 
1 — front bearing; 2 — central bearing, connection between two shaft parts 
(turbine shaft, compressor shaft); 3 — rear bearing; 4 — turbine; 
5 — compressor; 6 — central bearing journal 


a 


Fig. 2. Central bearing journal, connection between two shaft parts 
(turbine shaft, compressor shaft) 


The front bearing has a different number of rolling elements 
that the central bearing — that is why their characteristic sets 
are easily recognisable using the FDM-A method, Figs. 1 
and 2. The central bearing structure is the most heavily loaded 
subsystem of all three engine supports, due to most difficult 
cooling conditions and the highest deflection of the rotor unit. 
The outer ring of the central bearing is pressed into an elastic 
sleeve, the splines of which are fixed to the engine body. 


With unknown relations between the defect of the unit and changes in the depth of the generator output voltage frequency modulation. 
In fact, the presentation fi = f(t) has a discrete form as a set of points. In diagnostic practice, however, presentations are in use in which 


points are distributed with such a fine resolution that they give an impression of a continuous curve. 
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Occasionally, traces of rotation (slips) of the bearing ring in 
the sleeve can be observed during engine repairs. Also cases are 
recorded of cracks of the sleeve that fixes the central bearing 
outer raceway to the engine body. 

The inner bearing ring is pressed onto the compressor 
journal. The turbine shaft is introduced inside the journal by 
a splined connection. This connection makes it possible to 
transmit the torque from the turbine shaft to the compressor 
shaft. Longitudinal connection between these two parts is 
executed, in turn, by two half-rings and a special nut screwed 
onto the end of the central bearing journal. Breaking this 
connection may result in a disaster, as it disconnects the 
power turbine — there is no drive for pumps and other units 
necessary for proper engine operation. In this structure skews 
are frequently recorded which are the sources of additional 
dynamic loads, along with the forces generated by bearing 
assembly errors and radial clearances. 

From the point of view of journal load, very dangerous are 
large skew angles accompanied by small radial clearances in 
the central support — a source of heavy bending loads of the 
journal. The most dangerous cases of bearing operation include 
slips of the inner bearing ring on the journal, as they are the 
sources of not only of excessive dynamic forces, but also the 
increased level of temperature field, which can provoke journal 
plasticisation and reduce its mechanical strength. 

Due to the specific construction of the engine, of high 
importance it securing the coaxiality of the bearing seats during 
their assembly. Theoretically, all three bearings are to be situated 
along one axis, with a small permissible tolerance. Coaxiality 
changes lead to disturbances in bearing operation. Of high 
importance for each bearing is the state of the cage separating 
the rolling elements in the bearing. Its smooth motion guaranties 
correct motion of the rolling elements. The investigations 
were performed in two paths, Fig. 3, using signals from a DC 
generator (DC channel) and a three-phase AC rate generator (AC 
channel). These two paths complemented each other. 

The three-phase AC rate generator was a source of 
information on situations connected with slow-speed processes, 
including: 

% technical state of the power transmission box 

“* unbalance of the main shaft (with compressor and turbine), 
and ranges of radial clearances in bearing supports 

% skew of the main shaft — assessment could be made whether 
the skew of the compressor shaft with respect to the turbine 
shaft is increasing. As a final result, an assessment can be 
made whether this skew does not exceed a permissible limit 
defined by the bearing producer 

* crack of the sleeve that fixes the central bearing 

% increased ovalisation (circular shape error) of the bearing 
seat 

** number of cracked separation sections in the cage 

* error of journal perpendicularity with respect to the 
compressor disc rotation plane 

% increased torque on the engine main shaft. 


The DC generator provided the information about defects 
connected with high-speed processes. The processes which 
could be observed and diagnosed include: 
€ cooperation of the main shaft with the transverse vibration 

damper (central support) 
€ increased friction force in the front bearing as a result 

of longitudinal movements of the compressor shaft 

(h>N,p, 20.4) 


Measuring computer with FDM-A and FAM-C codes 


t 
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Fig. 3. Recording diagnostic signal from the examined engine 
using FDM-A and FAM-C methods 


€ locking of the front bearing cage (separator) — the appearance 
of the N.-th harmonic of the main shaft rotational speed 
(h= N) where N= number of rolling elements in the 
bearing 

© locking of the central bearing cage (separator) — the 
appearance of the N -th harmonic of the main shaft 
rotational speed (h = N) 

€ general level of wear of the bearing structure. 


PROBLEMS RELATING TO BEARING 
WEAR AND POSSIBILITIES OF THEIR 
OBSERVATION WITH THE AID OF FAM-C 
AND FDM-A METHODS 


A feature characterising the designs ofall types of bearings 
is the type of contact between the rolling elements and the ring 
raceways under load. This contact is of Hertzian type: linear in 
the roller bearings and of point nature in the ball bearings. It is 
decisive for determining the permissible load of the bearing® 
defined by its load capacity. The maximum load of the rolling 
elements resulting from the permissible Hertzian pressure in 
the rolling contact is a factor decisive for potential use of the 
bearing in a given case of load. According to Ref. [5] the highest 
probability that the machine defect will be detected is secured 
by the observation of the bearing structure, and this probability 
is equal to 70%. In laboratory conditions, the main cause of 
bearing defects was the fatigue wear of the working surfaces of 
their elements. According to literature data, the observations of 
defects of rolling bearings in real operating conditions indicate 
that the fatigue raceway defects do not exceed 36%, in general, 
of all causes of their withdrawal from operation, while in the 


According to ISO-76-1978 standard [6], permissible loads are the loads at which permanent deformation in the contact of the most heavily 


loaded rolling element with the raceway does not exceed 0,0001 of the rolling element diameter. In the USA, the standard ANSI/AFBMA 
std 9,11-1978 defines the maximum permissible limit of surface pressure as equal to 4000 Mpa. 
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remaining cases 55% of bearings lose their workability due to 
excessive friction wear, and further 9% - due to the corrosion 
and other destructive processes. 


Durability vs. load 


Fatigue durability of bearings depends on their load. As 
a consequence, observations of all instantaneous bearing 
overloads are crucial for the assessment of their lifetime. The 
FDM-A method makes it possible to trace, in an easy way, 
(high-speed) time-histories of the instantaneous velocity (DC 
channel), by observing the frequency, amplitude, and dynamics 
of retardations of the rolling elements. Indirectly, after gaining 
additional information from the slow-speed data recorded in 
the AC channel, we can also conclude about stresses taking 
place during rolling down the raceway. 

Hence, a conclusion can be drawn that the complete aircraft 
engines which, diagnosed using FDM-A and FAM-C methods, 
reveal lower (worsened) parameters cannot be mounted on 
airplanes designed for tasks connected with extremely high 
overloads, observed in acrobatic airplanes for instance. 


Friction wear 


It results from the investigations carried out by the research 
team headed by the author that changes of clearance dimensions 
can be monitored in the form of a number of parameters 
recorded using the FDM-A and FAM-C methods. These 
parameters include: 

æ» height of the characteristic set of the AC generator main 
shaft rotational speed first harmonic (AC measuring path) 

— total radial clearances 
~ value of the rolling coefficient p, defined as the ratio 

between the average angular velocity of the main symmetry 

axis of the bearing rolling element (roller) and the velocity 
of the main shaft. It is obtained from the quotient of the total 
frequency modulation band of the DC pulsation component 
in the bearing (DC measuring path) to the product of the 
rotational speed and the number of rolling elements in the 
bearing. The obtained p, value facilitates assessing the 
smoothness of the cage motion and makes it possible to 
determine the slip coefficient between the rolling elements 
and the raceway 

æ relative value of the AC path total pulsation — determining 
circumferential clearances between the cage window and 
the rolling element. 


Excessive bearing clearances are the sources of the raceway 
resonance vibration and provoke the instability of motion in 
the cage/rolling element system. The resonance vibrations have 
the high-frequency nature (of an order of a number of kHz) 
and can only appear when the SBK’ spectrum has the broad- 
band nature. The presence of slips in the rolling bearing 
differentiates the kinematics of particular rolling elements, 
which, as a result, leads to fluctuations of characteristic 
frequencies of the kinematic excitations. Due to the presence 
of clearances in the kinematic pair, each disturbance in the 
motion of the rolling elements is a source of additional impact 
forces in the cage/rolling element system. The volume of the 
energy introduced and lost during the collision depends on 
friction conditions between the contacting pairs composed of 


7 SBK - total shape error. 
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as the broad-band acoustic noise. 


main shaft speed. 


the rolling elements/cage partitions and the rolling elements/ 
bearing raceways. Intensive abrasive wear is accompanied, in 
turn, by higher bearing temperature, increased volume of wear 
products in the lubricating medium, and intensive ageing of 
this medium. Also the load distribution into particular rolling 
elements changes unfavourably, towards more heavily loading 
of the element which carries the highest force. 


Ageing phenomena 


Ageing phenomena are, perhaps, the least recognised 
issues. They require extremely long time intervals of engine 
observation, and precise records of their standstills and storage. 
It is well known that the fastest process of bearing ageing 
is observed on engines which are stored in unfavourable 
maintenance conditions, in high relative humidity (over n = 60%) 
and at large temperature changes. Especially harmful for 
the bearings are electrolytic corrosion processes, due to the 
use of various metals, with different alloy compositions 
(for instance, steel with different contents of carbon) and 
different electrochemical potential. It results from the author’s 
experience that the engines had smaller clearances after a few 
months of storage. 

The literature overview suggests that possible intensification 
of the ageing process depends strongly on changes in the 
operating temperature level caused by, for instance, high 
ambient temperature, and/or less effective cooling of the bearing 
due to insufficient mass flow rate of the cooling air or oil. 

A very dangerous type of ageing wear is also the 
abovementioned plasticisation of the bearing journal. 


SUCCESSIVE TYPICAL STAGES OF 
BEARING WEAR OBSERVED USING 
FAM-C AND FDM-A METHODS DURING 
THE OBSERVATION OF LONG-LASTING 
WEAR PROCESSES 


According to the practical diagnostic experience, successive 
stages (phases) of bearing wear [5] include noise wear, vibration 
wear, and thermal wear. 

In the examined engine the bearing cage is covered with 
a galvanic silver layer. In the first stage of bearing operation, 
i.e. the noise phase’, first we observe the decrease of the 
rolling resistance and the increase of radial clearances, due 
to wearing-in of the bearing. Then the silver is spalled out, as 
a result of which more and more heavy bearing operation is 
observed. The height of the characteristic sets of the bearings 
increases systematically, Fig. 6. The turning coefficient of the 
central rolling bearing’, defined as the ratio of the total band 
to the frequency of the N-th harmonic (where N — number of 
rolling elements), also increases in this time — the bearing cage 
rotates with increasing resistance. A recommended action in 
the noise phase is to observe and possibly correct the assembly 
errors in the engine main shaft/bearing structure system. In the 
above way such defects as the lack of journal perpendicularity 
to the rotation plane were detected in the investigations, a defect 
which manifested itself by increased depth of the amplitude 
modulation in the instantaneous frequency time-histories 
recorded in the AC generator path. In author’s opinion, as 
early as in this phase of bearing wear two different tendencies 
of turbine engine rotor bearing structure wear can be recorded, 


According to the literature [5] the dominating effects in the noise phase are macro-unevennesses of bearing elements manifesting themselves 


Rolling bearing turning coefficient — the ratio of the mean angular speed of the bearing rolling element (shaft) main symmetry axis to the 
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with two different presentation relations observed using FAM-C 

and FDM-A methods: 

a) type of structural-and-skew wear, in which systematic 
increase of the first subharmonic is observed in the 
AC generator path (increased skew) at relatively small 
increment of radial clearances in the central bearing. Most 
likely, at this type of wear large bending forces act on the 
journal, which in unfavourable situations can lead to its 
crack. The author did not observe such a phenomenon, 
and was only acquainted with the results of mechanical 
measurements done after the airplane crash. On the other 
hand, in the engines observed personally by him the bending 
forces did not lead to journal crack. They only provoked 
the increased load in the front bearing, with further multi- 
element resonance in the engine supports and strong 
destruction of the central bearing, Fig. 5. Certainly, these 
destruction phenomena refer to two next wear phases 

b) type of clearance wear, in which systematic increase of 
the first harmonic of the AC generator path is observed 
as aresult of rapid increase of the total radial clearances 
in the central support. At this wear form, the failure-free 
operation of the bearing is observed much longer than for 
wear type a. In some cases, an increased level of mechanical 
vibrations is observed for increasing time of operation. 
But the verifications which were carried out in that time 
did not reveal the state which would threat with bearing 
break-down. However, a literature overview suggests 
[16] that excessive clearances in rolling bearings can also 
be a source of excessive heat and lead to the destruction 
of the bearing structure. Indeed, an engine with bearings 
revealing excessive clearances has an unstable characteristic 
of changes of the load angle distribution into rolling 
elements at arbitrary manoeuvres of the airplane and that 
is why this state can be extremely dangerous, as intensive 
and uneven abrasive wear of the rolling elements can 
appear, accompanied by rapid increase of friction forces. 
In the areas of contacts between micro-unevennesses the 
temperature increases rapidly, reaching very high levels 
(flash temperature). As a result, the surface temperature 
of the contacting elements also increases, and so does 
their volumetric temperature. The temperature increase, 
especially in the micro-areas of friction contact, is a source 
of many changes, of both physical and chemical nature, in 
the materials used for bearing element and the lubricant. 


In the second stage, the vibration phase!®, systematic 
decrease of the height of the characteristic sets is observed due 
to the termination of the process of spalling out of the galvanic 
layer and washing out of the older spalling-out products, done 
by the lubricating oil system. As a result, the clearances between 
the bearing structure elements increase and the bearing motion 
resistance decreases. Also the rolling coefficient, defined in the 
literature [6] as the angle of the load distribution into particular 
rolling elements, decreases. At this stage, large clearances are 
observed in the bearing. The wear process is still in progress, 
but the mechanism of wear is controlled by the dynamic forces 
acting on the moving bearing elements in the environment of 
increased clearances. 

During the investigations a value of the rolling coefficient 
was recorded which does not make sense from the mathematical 
point of view. This means that not all rolling elements of the 
bearing roll down the raceway. In the vibration wear stage, 
situations can happen when the walls of the sleeve that fixes 
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observed in the literature on vibroacoustics [5]. 
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Fig. 4. Characteristic sets obtained from the DC generator path: 
1 — two characteristic subsets of the central bearing and the rear bearing; 
2 — characteristic set of the eighth harmonic; 
3 — characteristic set of the front bearing 
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Fig. 5. Characteristic sets obtained from the DC generator path 
for the state of general resonance in bearing structures 


the outer raceway to the engine body become thinner. This 
happens due to transient increase of forces locking the bearing, 
as a result of which the rotation-slip of the central bearing 
body (outer raceway) is observed in the sleeve, along with 
the erasing of the sleeve material due to sliding friction. The 
decreased thickness of the wall is observed in the FAM-C 
method as decreased amplitudes of the first harmonic and the 
first subharmonic. Most often this phenomenon is preceded by 
fading pulses, observed a large number of hours earlier using 
the FDM-A method in the DC generator path. Sometimes, 
due to remarkable reduction of thickness of the sleeve, or 
the appearance of large impact forces in the central support 
bearing, the sleeve cracks. Such a crack, Fig. 6, was detected 
using the FAM-C method in the three-phase AC rate generator 
path. In that time the height of the main shaft first harmonic 
characteristic sets decreased rapidly, by as much as nearly ten 
times, for all rotational speeds, except the range close to the 
maximum speed: n = 15100 rev/min (Fig. 7). 

It is well known [6] that the bearing structure becomes 
unstable after excessive deformations of bearing elements, 
shafts, and casing in operating conditions, or as a result of 
resonance vibrations, the reason for which is the coincidence of 
the exciting force frequency with the free-vibration frequency 
of the bearing. 

A very important parameter to be traced in the vibration 
wear Stage is the quality coefficient (Q) of the turbine engine 
rotor/bearing system [5]. The higher value of the quality Q, the 
smaller the damping in the system and the longer the transient 
processes. This way in case of large quality coefficients, each 
change of the nominal speed of the rotor, change of external 
dynamic forces acting on the rotor during manoeuvres of the 
airplane, and/or change of the load introduced by the units 
driven from the power transmission box, are the sources of 


Of highest importance in the vibrating phase are increased radial clearances, which result in increased amplitudes of mechanical vibrations 
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Fig. 6. Front bearing characteristic set height diagram for three successive 
time intervals of observation of the rolling bearing in the noise wear stage 
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Fig. 7. AC first harmonic characteristic set height changes during the crack 
of the sleeve that fixes the inner bearing raceway to the engine body: 


1 — after sleeve crack, 2 — before sleeve crack 
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strong nonlinear vibrations [15], accompanied by rubbing of 
the rotor elements against the casing and bearing seals, various 
strokes, etc. [14]. Determining the quality of the system is 
therefore of high importance for detecting possible causes of 
future failures. 

According to the literature [5, 14], the permissible quality 
for rotating machines is equal to Q = 5. With the aid of the 
FAM-C and FDM-A methods, determining Q is relatively 
easy as these methods present each bearing in the form of 
a separate characteristic set having a well defined amplitude 
and band width. During the investigations, a number of cases 
were recorded when some symptoms of mechanical resonance 
were recorded — when a characteristic set, or a number of sets 
decrease their frequency bands, and sometimes also increase 
their height (so called mechanical quality of the bearing 
structure increases, see Fig. 4 — detail 2). Sometimes these states 
refer to single elements of the turbine engine rotor centre, and 
manifest themselves as a “jump up” of a single characteristic 
set obtained using the FDM-A method. In other cases almost 
all characteristic sets increase their heights and decrease band 
widths. In those circumstances a multi-element resonance takes 
place and, according to [18], the system is in full resonance, 
Fig. 5. According to author’s experience, those states indicate 
advanced level of destruction of the bearing structure, Fig. 8. 
Before disassembly, this bearing revealed Q>11. During the 
disassembly of the bearing working in the resonance state for 
a relatively long time we can observe longitudinal indents on 
the raceways, like those reported in the literature. 

When the slip of the rolling elements over the bearing 
raceway becomes more intensive, the motion of the majority 
of the rolling elements may stop completely — the phase of 
the thermal wear begins. In this time the volumes of the heat 
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generated in the bearing structure are so large that they cannot 
be fully moved away by the oil system. As a result, the strength 
of the material structure of particular elements of the bearing 
decreases, which leads to excessive increase of windows in 
the separator cage, or even breaking the separation, Fig. 8. 
This state of bearing structure is presented in the instantaneous 
frequency time-histories f. = f(t) recorded in the AC path as 
repeatable, with respect to the shape and parameters, groups 
of pulses with the pulse numbers in a series proportional to 
the number of broken cells. Simultaneously, the pulsation 
amplitude increases from the level of about 10% up to 
160%. The increased temperature is also a source of thermal 
deformation of the central bearing ring, which easily looses 
clamp on the journal and is subject to slips. 


Fig. 8. Damaged central bearing, earlier diagnosed using the FAM-C 
method: 1 — broken cage separator ring 


Using the FDM-A method in the DC generator path we 
can precisely determine both the angle of these slips and the 
frequency of their occurrence. For this purpose the times of 
duration and the times of repetition of the fading pulses in the 
instantaneous frequency time-history f, = f(t) are observed and 
calculated. 

In author’s opinion, the most dangerous course of thermal 
wear stage, with respect to possible final effects, is that taking 
place at relatively small radial clearances!! in the central 
support bearing structure. In those circumstances, at large forces 
which tend to skew the two main shaft parts, the load is not 
effectively passed to the front support bearing and large part 
of the bending force is carried by the central bearing journal. 
All this can lead to the crack of the journal. 


ASSEMBLY ERRORS OBSERVED USING 
FAM-C AND FDM-A METHODS 


Bearing support assembly errors should be made distinctive 
and be analysed independently of the wear processes. They 
can act towards accelerating or delaying the time of action of 
particular wear stages during the bearing operation. This group 
of errors includes: 

- radial clearances of individual bearings 

- excessive circumferential clearances 

- error of misalignment of three bearing supports 
- error of turbine or compressor unbalance 

- error of bearing seat ovalisation 

- error of skew of the two shaft parts 

- error of compressor journal perpendicularity. 


Excessive radial clearances manifest themselves in the DC 
generator path as the increased amplitude of the characteristic 
set (presentation AF = f(f,)) of the central bearing. The excessive 


Small value of the radial clearance can be easily determined using the FAM-C method by detecting a relatively small amplitude of the 


first harmonic characteristic set on the presentation obtained from the AC generator path, and from small value (below 0.4) of the turning 


coefficient of this bearing, calculated using the FDM-A method. 
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circumferential clearance, in turn, manifests itself in the AC 
current generator path in the form of specific undercuts observed 
in the time-history f, = f(t). In general, all characteristics of the 
rolling coefficient are inversely proportional to the rotational 
speed, which can be easily explained as the effect of the 
increasing hydrodynamic forces which decrease the rolling 
friction resistance. Long-lasting observations, complemented 
by mechanical measurements, have led to a conclusion that if an 
inflexion point in the form of a minimum appears on the rolling 
characteristic, Fig. 9, this means remarkable misalignment 
between the bearing supports. 
0.70 
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Fig. 9. Changes of central support bearing rolling coefficient: upper curve 
— for negative misalignment, lower curve — without misalignment error 


The turbine or compressor unbalance error manifests itself 
in the presentations of the AC current generator paths as the 
increased amplitude of the characteristic set (presentation 
AF = f(f,) of the main shaft rotational speed first harmonic), 
accompanied by the widening of this band. If the bearing 
structure ovalisation error takes place, this set breaks down 
into two vertical characteristic sets. Skewing of the connections 
linking the two shaft parts manifests itself as increased height 
of the first subharmonic characteristic sets. The compressor 
journal perpendicularity error can be assessed by measuring 
the amplitude modulation depth in presentations f, = f(t). 


CONCLUSIONS 


The article presents the applications of the diagnostic 
methods FAM-C and FDM-A, which base on the analysis 
of output voltage frequency modulation parameters, for 
diagnosing technical states of turbine engine bearing supports. 
The use of these methods makes it possible to detect relatively 
early a defect in the diagnosed subsystem, before it becomes 
dangerous for the safety of airplane motion, for instance. 
Various assembly errors can be detected this way, and wear 
effects taking place in bearing structures can be monitored. 
Numerous parameters of the bearing operation can be traced, 
including: changes of friction resistance, changes of radial 
clearances in the bearing, changes of circumferential clearances 
between the cage and the rolling element, changes of the rolling 
coefficient. Moreover, it is possible to observe the phenomena 
taking place in the rotor unit, including: the skew between the 
turbine shaft and the generator shaft, and/or perpendicularity 
errors between the journal symmetry axis and the compressor 
disc rotation plane. Many hours earlier than when using the 
“traditional” methods we can observe, without interfering into 
the system, numerous wear phenomena that lead to dangerous 
failures of turbine engine bearing structures, such as bearing 
journal plasticisation, and/or disintegration of the separation 
gaps in the cage or of the entire circuit. Systematic use of the 
FAM-C and FDM-A methods makes it possible to determine, 
in a precise way, a safe time instant when the operation of the 
diagnosed object is to be stopped. The presented methods are 
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of general use and can be successfully applied for diagnosing 
technical states of turbine engines used in watercraft. 
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Problems of the starting and operating 
of hydraulic components and systems in low 
ambient temperature 


Part II 
Determining the clearance between cooperating elements during 
the hydraulic components start-up in extremely low ambient 
temperatures on the grounds of experimental research 


Ryszard Jasinski, Ph. D. 
Gdansk University of Technology 


ABSTRACT 


Severe winters and sweltering summers which more and more often occur nowadays are the 
reason why machinery designers face many difficulties when designing devices which will be 
serviceable in extreme ambient conditions. Hence, defining the principles and conditions of 
safe operation of hydraulically driven machines and devices is essential for their designers 
and operators. For this reasons the author did a series of tests of hydraulic component 
and systems in thermal shock conditions (cooled-down components were supplied with hot 
working medium). The experimental tests were carried out in the laboratory of the Chair of 


Hydraulics and Pneumatics, Gdansk University of Technology. They gave answer to the question how the 

effective clearance changes in such conditions and what parameters it depends on. Based on temperature 

graphs acquired from tests of heating up elements of hydraulic components it is possible to determine 
precisely change of clearance between cooperating elements. 


Keywords: hydraulic machines, hydraulic drives, diagnostics, hydraulic systems 


INTRODUCTION 


During the hydraulic components start-up in thermal 
shock conditions[ 1-9, 12-13] dynamic change of component’s 
elements dimensions due to thermal expansion is often 
encountered. Changes in the element’s clearings result in 
changes in the effective clearing between cooperating elements. 
Change of clearance is influenced, inter alia, by the temperature 
difference between hot oil and cold hydraulic unit, flow of 
working fluid, construction of the component, kind of material 
used and shape of component’s elements. 

Between the hydraulic component’s cooperating elements 
small clearances appear, measuring from few up to dozens 
um. Experimental research of various hydraulic components 
operating in low ambient temperatures undertaken by the 
author, allows to answer the question how the effective 
clearance changes in such conditions and what parameters does 
it depend on. The research included three groups of hydraulic 
components’ systems: supply system, control system and 
executive system. Fig. 1 presents the block diagram of typical 
hydraulic system exposed to sub-zero ambient temperatures. 
Main component of the supply system is the hydraulic pump, 
in the control system — valves and directional spool valve 
while in the executive system — hydraulic actuator (motor or 
cylinder). 


Using own methodology [1-8], the Author has tested 
many components and entire hydraulic system in thermal 
shock conditions. Tested components were placed in the low 
temperatures chamber in which -35°C temperature was reached. 
The components were supplied with the hot Total Azolla 46 
mineral oil in the temperature ranging from 20°C to 55°C. In 
the fixed and moving elements of components temperature 
sensors were placed, mainly thermocouples. The temperature, 
pressure and flow volume in the inlet and outlet channels of 
the components were also measured. 


HYDRAULIC PUMPS RESEARCH 


The object of the research were pumps of following 
various types: multi-piston axial pump of constant output with 
camshaft, multi-piston axial pump of variable output and two 
geared pumps. 

Researched camshaft pump of PWK 27 type (Fig. 2) of unit 
capacity 27.5 cm?/rev., operates at nominal pressure of 45 MPa, 
nominal speed of 1500 rpm and maximum 3000 rpm [19, 23, 
24]. The pump is characterized by high pressure-mechanical 
efficiency reaching 98 + 99 %, as well as by high overall 
efficiency at the level of 92 + 94 % [11]. 

Variable output multi-piston axial pump of PVM016R1D 
type (Fig. 3) is the smallest in the PVplus product line [21] and 
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Influence of the sub-zero temperatures 


Combustion 
engine or 
electrical 

motor 


Filters, 
coolers, etc. 


Fig. 1. Block diagram of machine 5 hydraulic drive with indication of the sub-zero temperatures influence on the operation of the components systems 


its characteristic parameters include: design output 16 cm?/rev., 
nominal pressure 28 MPa, max. pressure 35 MPa and rotation 
speed range between 300 and 3000 rpm. 


Fig. 2. Camshaft multi-piston pumps PWK-27 and PWK-78 
from Hydrotor company[10, 11,19] 


The arrangement consisting of PZ2-K-10 gear motor 
(Fig. 4) driving the PZ2-K-6.3 gear pump was tested. 

Determining the changes of effective clearance between 
cooperating elements of multi-piston pump in thermal shock 
conditions. 

In multi-piston radial and axial pumps or motors there are 
several places where loss of clearance between cooperating 
elements may occur. One of them is the effective clearance (1) 
between pistons moving in the cylinder block holes (Fig. 5). 
Change of clearance may be determined during start-up (t) 
according to formula: 


1,(t) = 1, + AL (x) -Al (1) 
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Fig. 3. Past-plate axial pump from Parker [21] company 


Fig. 5. Graphical interpretation of changes in the effective clearance 


It describes the changes of clearance height: caused by 
pressure — related deformation of cylinder block and piston 
Al, and by thermal expansion of both elements Al. 

Change of clearance height by the Al, value due to difference 
in linear thermal expansion of cooperating pump’s (motor’s) 
elements, in the effect of non-homogenous heating during 
start-up in thermal shock conditions may be described by the 
formula: 


Al(t) =h, : Br: [T,) -T,] + 


(2) 
z h, i Pa f (1) = 
where 
B,, 8, — linear thermal expansion coefficient’s of piston (T) 
and cylinder block (B) 
T,(t) — piston temperature 


T,(t) — cylinder block temperature 

— measuring temperature of part’s linear dimensions 

— linear diameter dimension of: piston (T) and 
distribution bush openings (B). 


PREPARING THE PWK27 PUMP FOR 
MEASURING THE CHANGE 
OF ELEMENTS’ TEMPERATURE 


In the tested pump PWK (Fig. 2), built of two groups of 
elements: fixed (Fig. 6) and mobile (Fig. 7) the thermocouples 
were placed. 

Mobile elements of the pump include: piston, screw 
throttle, hydrostatic slipper, distribution bush, resistance shield, 
clamping ring, separator, shaft, cam. 

Fixed components of the pump include: cover, front 
housing, cylinder block, cylinder block cover, rear housing, 
guiding bush, socket. 

Constant of temperature measurements allowed 
determination of the clearance between cooperating elements 
(pistons and cylinder bushes) of the pump during start-up in 
thermal shock conditions. 

For the sake of researching the process of piston and 
slipper heating, they were purposely made with bores for fixing 
thermocouples. (Fig. 7, 8). Thermocouple leads placed in the 
slipper and piston were lead through purposefully made bores 
in the separator and cylinder block (Fig. 9). 


Fig. 6. Distribution of thermocouples (T1-T6) in the holes bored in PWK27 
pump 5 fixed elements: T1 — cover, T2 — front housing, T3 — cylinder block 
— area no. 1, T4 — cylinder block — area no. 2, T5 — cylinder block cover, 
T6 — rear housing 


Fig. 7. Thermocouples in PWK27 pumps mobile parts: 
a) piston, b) hydrostatic foot 


Fig. 8. Piston unit adapted to measuring the heating process 


POLISH MARITIME RESEARCH, No 1/2009 63 


Fig. 9. Connection of thermocouples leads to the piston 
and hydrostatic foot in researched PWK27 pump 


RESEARCH OF PWK27 PUMP IN 
THERMAL SHOCK CONDITIONS 


Experimental research of PWK27 pump operation during 
start-up in thermal shock conditions were executed according 
to methodology [1] designated for the following start-up 
conditions: 

* oil temperature: ~ +48°C 

* ambient temperature: -21 + +23°C 

* pump’s shaft rotation speed: 500 + 2500 rpm 
* pump’s load: 4.5 + 12 MPa. 


Several measurement series of PWK27 pump start-up 
in thermal shock conditions were run. The results include, 
inter alia, characteristics of temperature change of: heating 
fixed elements, temperature of incoming and outgoing oil, 
temperature of leakages and characteristics of torque, oil 
pressure on the outlet, output, external leakage volume, 
rotation speed of the shaft, volumetric, torque and overall 
efficiency. 

One start-up test of PWK27 pump was conducted for 
the following parameters: outlet pressure 9.15 MPa, output 
30 dm?/min, pump’s shaft rotation speed 1100 rpm, oil 
temperature 49°C, ambient temperature -12°C. 

On the grounds of recorded temperature graphs (Fig. 10) it 
is possible to determine that the fastest heating element in the 
pump is the cylinder block unit consisting of cylinder block 
and its cover. 
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Fig. 10. Temperatures in pump 5 fixed parts (Fig. 6) 
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The fastest temperature build-up in the cylinder block occurs 
near the pump’s internal canals that carry the main stream of hot 
oil and near the pump’s camshaft (thermocouple measurement 
no. 4). Almost as quickly the temperature rises in the common 
areas of piston and guiding bush (thermocouple no. 5). Much 
slower however the temperature rises in front (thermocouple 
no. 2) and rear housing (thermocouple no. 6). In the first 200 
seconds of pump’s operation temperature graphs of the two are 
identical, only after that the rear housing begins to heat much 
quicker. It is this element that reaches the highest temperature of 
all the pump’s fixed elements under set conditions of operation. 
Whereas the slowest heating up element is the pump’s cover 
(thermocouple no. 1). 

It is worth considering that same of PWK pump’s fixed same 
elements (Fig. 10) while operating under constant pressure reach 
higher temperatures than this of the hot oil feeding the pump. This 
is exemplified by the comparison of characteristics on Fig. 11 
with those on Fig. 10. It is the effect of friction heat generation 
in the stream of fluid flowing through pump’s crankcase and 
mixed by the resistance shields of feet and pistons. 


= oil temperature on outlet 
= oil temperature on inlet 
A leakage temperature 


Temperature [°C] 
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Fig. 11. Oil temperatures: on pump ° inlet, 
on pump 8 outlet and external flows 


In the initial period of pump’s operation the temperature 
of leaking oil rises far slower than the temperature of main 
working stream on the inlet and outlet (Fig. 11). However, after 
some time the temperature of leakage reaches much higher 
temperature than the feeding stream temperature. 

During tests PWK27 pump reached very high volumetric 
efficiency at least 0.98 (Fig. 12). 


0 100 200 300 400 500 600 700 800 900 


Fig. 12. Overall (n ), volumetric (ņ ) and torque (n,,,) efficiency for the 
following parameters: working pressure 9.15 MPa, output 30 dm*/min, 
ambient temperature -12°C, oil temperature 49°C 


Torque efficiency remains at the level of 0.85. It largely 
depends on pump’s load that in the considered situation only 
amounted to 9 MPa. The higher is the working pressure of the 
pump, the higher its torque efficiency, and in the same time the 
lower its volumetric efficiency. 


DETERMINING THE EFFECTIVE 
CLEARANCE BETWEEN THE PISTON AND 
CYLINDER OF THE PWK 27 PUMP 


On the grounds of analysis of several heating cycles of the 
piston achieved in few series of tests it was determined that 
the rate of piston heating is always lower or equal to the speed 
of the hydrostatic slipper heating. Fig. 13 presents the heating 
processes of the piston and the cylinder block cover. 
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Fig. 13. Changes in temperature of piston hydrostatic slipper and cylinder 
block and difference of temperatures between piston and cylinder block for 
start-up parameters: ambient temperature -20°C, oil temperature 48°C, 
rotation speed 1100 rpm, pressure above piston 5.2 MPa 


Changes of temperature result in changes of radial clearance 
shown in Fig. 14 for pump cooled to -20°C supplied with hot 
oil +48°C at shaft rotation speed 1100 rpm. 
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Fig. 14. Change of clearance between the piston and distribution bush 
for the parameters of pump 5 operation after start-up: ambient 
temperature -20°C, oil temperature 48°C, rotation speed 1100 rpm, 
pushing pressure 5.2 MPa 


The most significant reduction of clearance occurred after 
ca. 30 seconds after pump’s start-up (Fig. 14). It is this period 
when the biggest difference in temperatures occur between 
the temperature of the piston and distribution bush. Further 
on, during start-up the temperatures of cooperating elements 
become even resulting in increase of radial clearance to desired 
value. 


RESEARCH OF THE GEAR PUMP PZ2-K10 
OPERATING AS MOTOR 


In the pump or gear motor the reduction of radial clearance 
between teeth points and internal surface of the housing may 
occur. 


The lack of axial clearance however will not occur because 
the pump includes mobile elements providing compensation for 
axial clearances. Also the jamming of bearing’s casing made of 
steel in the housing made of aluminium will not occur. Based 
on geometric measurements it was determined that the radial 
assembly clearance in the tested pump PZ2-K10 (Fig. 15) 
amounts to 0.087 mm in the ambient temperature of 20°C. 


Fig. 15. Places of temperature measurement in gear pump PZ2-K10 
elements: plate, cover, gear wheel, housing 


According to the equation | the effective clearance not 
only depends on assembly clearance but also on deformation 
of elements due to influence of pressure and different thermal 
expansion on cooperating elements. The higher the oil pressure, 
the bigger is the slack of gear shaft and its displacement inside 
the bearing. This in turn affects the dimensions of gaps between 
cooperating elements. It has been assessed that at the maximum 
operating pressure of the tested pump the displacement of gear 
shaft in the bearing and its slack may reach 21 um [17, 19]. This 
would cause the decrease in dimension of gap between gear 
teeth tips and hole in the housing on the suction space side. 

The change in element’s dimensions due to difference in 
thermal expansion of gear pump elements with housing made 
of aluminium alloy in thermal shock conditions is described 
by the following equation: 


Al(t) = 11+ 10%- h, - [T,(t) — 293] + 
-22.3 - 10% - h, + [T,(t) — 293] [m] 


— linear dimensions: tip diameter of gear tooth (Z) 
and housing opening diameter (K) [m] 

T,(t),T,(t) — temperature of gear wheel and housing in the 
process of heating up [K] 

— temperature of linear measurement of components 
[K]. 


293 


Since the linear thermal expansion coefficient for aluminium 
alloy is almost twice that of steel when lowering the ambient 
temperature the initial clearance will also decrease. (Fig. 16). 
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Fig. 16. Initial radial clearance in PZ2-K10 gear pump (motor) 
depending on ambient temperature 


Based on experimental research, the courses of temperature 
changes in housing and gear wheel of PZ2-K 10 pump (operating 
as motor) were determined for the flow volume in the range 
of: from 4.5 dm*/min to 14 dm?/min and motor’s load: from 
0.5 Nm to 17 Nm. 

One of the tests was conducted for the following parameters: 
oil temperature 48°C, ambient temperature -16.8°C, working 
pressure 5 MPa, pumps rotation speed 1020 rpm. 

It is to be noticed, that the gear wheel heated up quicker 
than the housing. Maximum difference in temperatures between 
heating elements of the unit reached 12°C in the 80th second 
after pump’s start-up (Fig. 17). 
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Fig. 17. Temperature of housing, gear wheel and difference in temperatures 
between those elements of PZ2-K10 pump operating as motor. 


Based on temperature changes graph acc. to Fig. 17 it is 
clear that during gear pump start-up the rapid decrease of radial 
clearance between teeth tips and internal surface of housing 
does not occur (Fig. 18). In case of hydraulic components 
made of aluminium alloys (housing) and steel (wheel) the 
disappearance of radial clearance doesn’t occur. 
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Fig. 18. Change of effective clearance between tips of teeth and internal 
housing of PZ2-K10 pump operating as motor, during start-up under 
following conditions: oil temperature 48°C, ambient temperature -16.8°C, 

working pressure 5 MPa 
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In case of pumps and gear motors with housings made of 
cast iron the changes of clearance will look differently than 
on Fig. 18. Such components will be susceptible to thermal 
shock conditions. 


RESEARCH OF INADEQUATE OPERATION 
OF DIRECTIONAL SPOOL VALVES IN 
THERMAL SHOCK CONDITIONS 


When in a hydraulic circuit a sudden case of hot oil 
flowing through cold valve (thermal shock conditions) 
occurs the uneven heating of cooperating elements of spool 
and casing will follow. The spool is much smaller than 
the casing, therefore it’s thermal capacity is far smaller. In 
comparison to the casing it has relatively large surface flushed 
by passing oil. Therefore it will heat up much faster. With 
rising temperature of the spool its dimension increase faster 
against the opening in the casing leading to the decrease 
of the clearance. The clearances are already quite small in 
order to minimize leakages. Too small clearances between 
the spool and casing may result in substantial friction forces 
hindering or even preventing their movement against each 
other. Depending on the construction of the component and 
its maker the clearances between the spool and casing fit in 
the range of 5 - 25 um [18]. 

Effective clearance (1,) between casing and spool (Fig. 19, 
20) depends on assembly clearance (1,), elastic distortions 
of spool valve’s elements due to pressure (Al ), difference in 
thermal expansion of various components (A1). 


Start of supply 


Fig. 19. Change of effective clearance between spool and casing 
of cooled spool valve fed with hot oil during start-up 


In thermal shock conditions the operation of the spool valve 
is seriously affected by the thermal expansion characteristic 
of materials used for making specific elements of the spool 
valve. Spools of hydraulic spool valves are made of steel while 
their casings are usually of cast iron. The linear expansion 
coefficient for cast iron equals 10.5*10° 1/K, while for steel it is 
11*10° 1/K [14, 15]. 

Taking under account the influence of thermal expansion 
coefficients the effective clearance may be calculated from the 
relation of (1) and (2). 

Based on numerical calculations it was assessed that 
the influence of pressure on effective clearance is of lesser 
importance that the influence of substantial disproportion in 
temperature of cooperating elements (spool and casing). The 
confirmation of the fact is found through the calculations of the 
difference in thermal expansion Al, and distortions resulting 
from the influence of pressure Al performed for the spool valve 
RE2510/101 from Hydrotor (Fig. 21), where the spool diameter 
measured 18 mm. With the spool and casing temperature at 
AT=20°C the assembly clearance will decrease by Al, =3.9 um, 


for the temperature difference AT = 45°C it will decrease by 
Al, = 8.7um. The influence of oil pressure (16 MPa) on 
components of the spool valve will only cause the assembly 
clearance to change by Al, = 1.5 pm. 


Points of common P 
operations of spool i i 
and casing 


Fig. 20. Points of common operation of spool and casing 
and graphical interpretation of change of effective clearance 


SPOOL VALVES TESTED IN LOW 
AMBIENT TEMPERATURES 


In the lab of the Chair of Hydraulics and Pneumatics 
Gdańsk University of Technology following valves for control 
the pressure and flow direction were tested in thermal shock 
conditions [5, 6]: 

+ relief valve type UZPP16 from PONAR WADOWICE 
company 

relief valve type DBW 20 A2-5X/350Y S6EG24N9Ķ4R12 
from REXROTH company 

dual stage spool valve 4WEH16C33/6AW220-50 from 
REXROTH company 

electro hydraulically controlled spool valve RE2510/101 
from HYDROTOR company (Fig. 21) 

proportional spool valve PVG 32 from SAUER DANFOSS 
company (Fig. 22) 

servovalve 4WS2EM10 - 45/20B2T315Z8EM from 
REXROTH company (Fig. 23). 


+ + + + + 


Hydraulic units are often fitted with electro-hydraulically 
controlled spool valve. Therefore many constructions of 
electrically controlled hydraulic valves were tested. 


b, 


Fig. 21. Electro-hydraulic spool valve RE2510/101 [19] 


Fig. 23. Servovalve 4WS2EM10 - 45/20B2T315Z8EM [18] 


Tested sliding spool valve belong to two groups depending 
on their function: 


I — spool valve used for controlling the flow direction - 
classical (Fig. 21, 24) 
II — proportional spool valve and servovalves additionally 


used to control the volume of flow (Fig. 22, 23, 25, 26). 


Nominal diameters of spools in case of PVG 32 and 
RE2510/101 spool valve amounted to 18 mm. Initial 
clearance between spool and casing in PVG 32 spool valve 
amounted 7um, while in RE2510/101 spool valve it was 
6um (measurements performed at Metrology Lab of Gdansk 
University of Technology). 

In the body of RE2510/101 spool valve (Fig. 24) bores 
were made for placing thermocouples. Bores for thermocouples 
T7 and T8 are placed in axis but on different depths, near the 
source of flowing oil while the bore for thermocouple T3 is 
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away from the hot oil channel. In order to test the spool two 
bores were made to accommodate thermocouples. The bores 
were made in axis however to different depths. One of the 
bores was located near the place where cylindrical surface of 
the spool was flushed with hot oil. The other was located at 
the spool ’s front. 
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Fig. 24. Spool valve RE2510/101 from Hydrotor company T1, T2 — 
thermocouples placed in spool T3, T7, T8 — thermocouples placed in casing 


Fig. 25. Spool valve PVG 32 from Sauer Danfoss. T1, T2- thermocouples 
located in the spool T9, T10, T11, T12 thermocouples located in the fixed 
elements of the spool valve 
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Fig. 26. Servovalve 4WS2EM10 - 45/20 from Rexroth with temperature 
(thermocouples) and pressure sensors 
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For the purpose of measuring the heating up of PVG 32 
spool valve the bores were made in spool and casing in which 
thermocouples were inserted. Their layout is presented on Fig. 
25. In the spool a bore was made reaching half of its length. 
T1 thermocouple was inserted into the full depth of the bore 
while the T2 was placed in its beginning. T9, T10, T11 and 
T12 thermocouples were placed in fixed elements of the spool 
valve. 

The diameter of the researched servovalve 4WS2EM10 
is 6,495 mm, while the clearance between the spool and 
cylindrical bush in the casing is 5 um. 

In the fixed elements of 4WS2EM10 - 45/20 servovalve 11 
thermocouples were placed (Fig. 26), allowing for presentation 
of heating processes for various states of supply. Through 
pressure indication, the pressure sensors (Fig. 26) enabled the 
determination of reaction time to the control signal and change 
of this time with changing ambient temperature. 


SAMPLE RESULTS OF TESTING PVG 32 
SPOOL VALVE 


Several tests series of PVG 32 spool valve were conducted 
for the temperature differences between hot oil and cooled unit 
kept in the range of AT = 20+ 75 K. 

One test of the spool valve cooled to 0°C was conducted at 
abrupt changes of flowing oil pressure from zero to 32 dm?/min 
at temperature 52°C. The difference in temperature between 
the spool and casing was determined (Fig. 27). Maximum 
difference in temperature between the body and the spool during 
heating reaches 27.9°C and occurs in the 30th second. In this 
time the temperature of casing is 13.8°C and spool ’s 41.7°C. 
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Fig. 27. Average temperatures of casing and spool and the difference 
in temperatures between the two during feeding: oil temperature 52°C, 
ambient temperature 0°C, flow volume 32 dm*/min 


Using the formula (1) under assumption that the distortion 
of the spool and the casing due to pressure are negligible, the 
effective clearance between the components was determined 
(Fig. 28). 


Clearance [um] 


Time [s] 
Fig. 28. Effective clearance between external surface of the spool 
and internal surface of the casing during feeding: oil temperature 52°C, 
ambient temperature 0°C, flow volume 32 dm*/min 


Based on the effective clearance graph (Fig. 28) it is possible 
to determine that the disappearance of clearance between spool 
and casing due to difference in thermal expansion between the 
two will not occur. 


RESEARCH OF THE LOW SPEED 
HYDRAULIC MOTORS 


The first part of the article [8] deals with symptoms of 
correct and incorrect operation of hydraulic motors in thermal 
shock conditions based on pressure and turning speed. The 
results of tests of satellite hydraulic motor were presented. 

Also low speed orbital hydraulic motors without axial 
clearance compensation were tested in low ambient temperatures. 
Such motors are commonly used in driving arrangements of 
machines. 

Researched satellite motor SOK 100 from ZUO Hydroster 
[20] (Fig. 30) of geometric working capacity 0.4 dm?/rev. may 
operate at the range of speed from 0 to 160 rpm (n,a = 250 rpm) 
and nominal pressure p, = 16 MPa (P am 25 [MPa]). Assembly 
axial clearance was 23 um. 

Orbital motor GMR 160 from Rexroth (Fig. 33) [18] of 
geometric working capacity 162.9 cm*/rev according to its 
maker is characterized by nominal speed of 375 rpm (nax 
= 470 rpm), nominal pressure p = 12 MPa and max. pressure 
Pmax = 22-5 MPa. 

TF170MW260 motor from Parker company[21] of 
geometric working capacity 169 cm*/rev according to its 
maker may be operated continuously up to the nominal speed 
of 355 rpm and nominal difference in pressure (between in and 
outtake) 13.8 MPa, and momentarily at n „= 440 rpm and max. 
pressure difference 20.7 MPa. The maximum supply pressure 
may reach 30 MPa. 


PREPARATION AND SAMPLE RESULTS 
OF TESTING THE SOK 100 MOTOR 
(FIG. 29) CONDUCTED IN THERMAL 
SHOCK CONDITIONS 


The motor was made to specific order at ZUO HYDROSTER 
in Gdansk. In the fixed elements specials bores for thermoelectric 
sensors were made for measuring temperature in front and 
rear covers and the housing. Following that, the motor was 
assembled with semiconductor temperature sensors (AD590) 


placed in mobile elements (two sensors in the rotor and one 
in the shaft). 
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One test of SOK 100 when the motor operated incorrectly 
was conducted for the following conditions: flow volume 
100 dm?/min, initial temperature -20°C, oil temperature 50°C. 

Start-up of motor in this conditions caused substantial rise of 
temperature in mobile elements of the motor. The satellites heated 
quicker than the rotor. Fig. 30 presents the graphs of temperature 
in rotor and housing. The substantial rise in rotor’s temperature is 
clearly visible (Fig. 30). This, in turn results in great difference in 
temperatures of the rotor and housing leading to disappearance 
of clearance between the two (Fig. 31). Substantial rise in 
temperature caused the dry friction between the elements. 

In result the friction heat appeared directly effecting elements 
of me motor. In turn satellites and rotor heated up quicker. 
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Fig. 30. Temperature of rotor and housing with difference in temperatures 
between the two in SOK 100 against timeline 
25 
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Fig. 31. Change of effective clearance between the covers and rotor of SOK 
100 during start-up in thermal shock conditions 
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Fig. 29. Motor SOK: 1 — housing; 2 — rotor; 3 — satellite; 4 — shaft; 5, 6 — covers; 7 — setting ring; 
8, 9 — non-return valve; 10, 11 — rings; 12 — pin; 13 — manifold 
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Radial clearance 


Axial 


Fig. 32. Orbital motor GMR 160 with temperature sensors: 1 — shafting, 2 — front cover, 3 — housing, 4 — geared shaft, 5 — rotor, 6 — body, 7 — rear cover, 
8 — distribution disc, 9 — roller (markings TI + T8 indicate place of temperature measuring using thermocouples of the same marks) 


PREPARATION OF ORBITAL MOTORS 
FOR TESTING IN THERMAL SHOCK 
CONDITIONS 


The orbital GMR 160 motor (Fig. 32) includes construction 
nodes susceptible to disappearance of clearance in thermal 
shock conditions. 

The nodes are as follows: 

e axial gaps between surfaces of rotor, cover and distribution 
disc 

e gaps between perimeter surfaces of rotor’s and housing’s 
gear 

e gaps between external surface of camshaft and internal 
surface of housing. 


In order to assess specific elements of GMR 160 motor in 
various places 8 temperature sensors (thermocouples) were 
placed. In the rear cover thermocouples T7 and T8 were placed. 
Six thermocouples were inserted into the bores in the housing 
- thermocouples T1, T3, T5 (away from the heat source) and 
T2, T4, T6 (close to the heat source). 

The second orbital motor tested in thermal shock conditions 
was TF170 (Fig. 33) with frontal distribution. 

The motor was prepared for tests by making bores in fixed 
elements for thermocouples, i.e. in spool valve, housing near 
the working chamber, housing near the roller, body near the 
oil channel, body away from the heat source. Sensors were 
not placed in mobile elements due to complicated kinematics 
of their operation. 

Orbital motors GMR 160 are TF 170 are differentiated 
by the distribution and the flow of oil through the channels 
resulting, in turn, in different way of heating up motor’s 
elements. 


Fig. 33. Orbital motor TF 170: 1 — spool valve, 2 — housing, 
3 — geared shaft, 4 — roller, 5 — rotor, 6 — body 
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INFLUENCE OF ORBITAL MOTORS’ 
CONSTRUCTION ON THEIR OPERATION 
IN THERMAL SHOCK CONDITIONS 


In order to assure proper operation of hydraulic system in 
thermal shock conditions suitable hydraulic motors, resistant 
to such conditions have to be used. Such requirements are well 
met by the tested orbital motor TF 170 from Parker company 
(Fig. 33, 34). In its housing there are additional channels 
for flowing oil thus enhancing the surface of heat exchange 
resulting in quicker heating up of the housing. It was also 
found out that with such construction (Fig. 35) the maximum 
difference in temperatures between the housing and rotor didn’t 
exceed 20°C during start — up under thermal shock conditions 
for the temperature difference AT „ of oil and surrounding 
reaching 77°C. 
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Fig. 34. Components of GMR 160 orbital motor (on the left) 
and TF 170 (on the right). In the TF 170 motor additional 
channels for flowing oil are provided 
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Fig. 35. Difference in temperature between rotor and housing of TF170 
motor (volume of flow 64 dm?/min, initial temp. of motor -27°C, oil temp. 
50°C) and GMR160 motor (volume of flow 70 dm*/nin, 
initial temp. -26°C, oil temp. 52°C 


Orbital motor GMR 160 (Fig. 32 and 34) provides no 
additional heat exchange surface from oil to housing. It was 
also found out (Fig. 35) that the maximum difference in 
temperatures between the rotor and housing reached 46°C for 
AT , 78°C. 


Having the experimental data as per Fig. 35 and based on 


the relation (1) change of radial clearance was determined for 
motors GMR 160 and TF170 (Fig. 36). 
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Fig. 36. Graph of radial clearance in TF170 and GMR160 motors 
(start-up conditions as in case of Fig. 35) 
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Radial clearance in TF170 motor, between the rotor 
and housing decreased maximum by 12um, while in case 
of GMR160 motor, under similar operating conditions the 
clearance decreased by 271m and was close to vanishing 
completely (Fig. 36). 

Thanks to substantial initial clearance the GMR 160 motor 
proved resistant to thermal shock conditions and operated 
correctly at temperature difference AT , | of oil and surrounding 
up to 78°C. 

In case of TF170 motor from Parker during start-up under 
thermal shock conditions only slight change of radial clearance 
between the rotor and housing occurs (Fig. 36). This signifies 
substantial resistance of such motors during operation under 
given conditions. 


SUMMARY 


O The main construction issue decisive for admissible supply 
parameters of cooled hydraulic components in thermal 
shock conditions is the effective clearance. 


O Based on temperature graphs acquired from tests of 
heating up elements of hydraulic components it is possible 
to determine precisely change of clearance between 
cooperating elements. 


O There are many factors influencing change of effective 
clearance, i.e.: initial clearance, initial temperature of 
the component, temperature of feeding oil, volume of 
flow, material used for making the unit and element’s 
construction. 


O Using suitable materials for cooperating elements of 
hydraulic components it is possible to protect them fully 
against thermal shock conditions. This is well exemplified 
by high speed pumps with housing made of aluminium 
alloy and gear of steel. 


O During start-up in thermal shock conditions the supply and 
executive components present low torque efficiency and 
therefore low overall efficiency. This is due to substantial 
friction in working agent flow through the component’s 
channels. One example is the torque efficiency of PWK 
27 pump determined for start-up at temperatures -20°C 


and +24°C. During start-up in thermal shock conditions 
the pump presents decreased torque efficiency (Fig. 37). 
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Fig. 37. Hydraulic — mechanical efficiency of start-up of PWK27 pump 
Jor the following conditions: ambient temp. t,,, = +24°C and t „= -20°C, 


revolution speed n = 1100 [rpm], oil temp. t, = 48°C, 
working pressure 5 MPa 


O Apart from above described experimental method of 
determining parameters for correct start-up of hydraulic 
components in low ambient temperatures, Author uses and 
recommends two other methods: analytical and computer 
aided simulation. With the use of these methods one may 
determine, with suitable precision, parameters for correct 
start-up of hydraulic system in thermal shock conditions 
without the necessity of using costly, troublesome and 
often destructive experimental researches. These methods 
will be described in following components of this article, 
concerning issues addressed here. 
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Estimation of the probability of propulsion loss 
by a seagoing ship based on expert opinions 


Alfred Brandowski, Prof. 
Gdynia Maritime University 


ABSTRACT 


2 The event of the loss of propulsion function has been defined as hazardous event to 
a seagoing ship. It has been formalized. The procedure of acquisition of expert opinions 
on frequency of the event occurrence has been described. It may be considered to be of 
a numerical-fuzzy character. The fuzzy part was transferred to the numerical form by the 
pair comparison method. An example of the ship propulsion system comprising a low speed 
internal combustion engine and a fix pitch propeller illustrates the method presented. It may 
be used wherever a hazard analysis has to be performed of a system involving human and 


technical aspects and there is a shortage of objective data on the investigated object. 


Keywords: ship, loss of propulsion, probability, estimation, expert 


INTRODUCTION 


Loss of the propulsion function by a ship is one of the most 
serious categories of hazardous events! in shipping. In specific 
external conditions it may lead to a loss of ship together with 
people aboard. The loss of propulsive power may be an effect of 
the propulsion system (PS) failures or of errors committed by the 
crew in the system operation process. In the safety engineering 
language we say that the propulsion loss probability depends 
on the reliability of the PS and of its operators. Determination 
of that probability is in practice confronted with difficulties 
connected with shortage of data on that reliability. This 
pertains particularly to the cases of estimation in connection 
with decisions taken in the ship operation. In such cases we 
have to rely on subjective estimations made by persons with 
practical knowledge in the field of interest, i.e. experts. The 
experts, on the other hand, prefer to formulate their opinions 
in the linguistic categories, in other words in the language of 
fuzzy sets. The author’s experience tells also that in the expert 
investigations it is difficult to maintain proper correlation 
between the system data and the system component data. 
The paper presents a method of the subjective estimation of 
propulsion loss probability by a ship, based on the numerical- 
fuzzy expert judgements. The method is supposed to ensure 
that proper correlation. It is adjusted to the knowledge of 
experts from ships’ machinery crews and to their capability of 
expressing that knowledge. 


The method presented has been developed with an intention 
of using it in the decision taking procedures in risk prediction 
during the seagoing ship operation, in the shortage of objective 
reliability data situations. 


DEFINITION OF THE SHIP PROPULSION 
LOSS AS A HAZARDOUS EVENT 


The propulsion hazard is connected with the loss by the PS 
system of its capability of performing the assigned function, i.e. 
generating the driving force of a defined value and direction. It 
appears as an effect of a catastrophic failure” of the PS. Such 
failure may cause immediate (ICF) or delayed (DCF) stoppage 
of a ship. In the latter case the stoppage is connected with 
renewal, which may be carried out at any selected moment. 
It is obvious that only the former case of the forced stoppage 
creates a risk of damage or even loss of ship - it is a hazardous 
event. 

We will relate the probability of ICF to an arbitrary time 
interval determined by the analyst. For instance, it may be 
duration of one trip, time interval between the ship class 
renewal surveys or one year, as it is usually assumed in risk 
analyses. Such an approach is useful in the ship operation risk 
management process. 

The ICF type failure consequences may be divided into 
casualties and incidents (IMO 1997). In general, the ship 
casualties are non-repairable at sea by means of the ship own 


' Hazardous event is defined as an event bringing about damage to human beings as well as to the natural and/or technical environment. 


It is also called “accident” or “initiating event”. 


? Catastrophic failure is defined as loss of the capability of performing by the object of its assigned function. 
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resources and may have very serious consequences, with the 
ship towing at the best and the loss of ship at the worst. The 
problem of consequences is not the subject of this paper. 

The ICF type failure frequency depends mainly on the type 
of PS and the ship operation mode (liner trade, tramping etc.). 
On the other hand, the consequences are strongly dependent on 
the ship size and type and the environmental conditions, first 
of all the water region, season, time of day, atmospheric and 
sea conditions. They are also dependent on the navigational 
decisions and on the type and fastening of cargo in the holds 
and on deck. In general, these are the factors connected with 
the type of shipping carried out and the shipping routes the 
ship operates on. 


FORMAL MODEL OF ICF EVENT 


We assume the following: 

ve We are interested only in the "active" phase of ship 
operation, when it is in the shipping traffic. We shall 
exclude from the model the periods of stays in ship repair 
yards or in other places connected with renewals of the 
ship equipment. 

ve The investigated PS system may be only in the active usage 
or stand-by usage state. The ICF type PS failures may occur 
only in the former state. 

ve A formal model of the ICF type PS failures is the 
homogeneous Poisson process (HPP). This assumption is 
justified by the expert elicited data, which indicate that this 
type of failures occur fairly often, several times a year, but 
their consequences in general mean only a certain loss of 
operation time. More serious consequences, causing longer 
breaks in the normal PS system operation, occur seldom. 
The exponential distribution of time between failures, taken 
place in the HPP stream model, is characteristic of a normal 
operation of many system classes, including also the ship 
systems [1, 4]. It is appropriate in the case when the modeled 
object failures and the operator errors are fully random 
abrupt failures and not gradual failures caused by the ageing 
processes and/or wear of elements. This corresponds with 
the situation when scrupulously performed inspections and 
renewals prevent the latter type of failure from occurring. 

ve Experts are asked only about two numerical values: 
number of ICF type failures N(t) during time period 
t= 1 year (8760 hours), and the time at sea percentage share 
«100% during their seamanship period - this is within their 
capability of answering. 

ve The opinions on the failures of PS system components are 
elicited in the linguistic form. 


The seagoing ship system active usage time t® is strongly 
correlated with the specific ship operational state times, mainly 
with the “at sea” state including “sailing”, “manoeuvres” and 
“anchoring”. The following approximation may be adopted for 


the system, also for the PS: 


t® =t™ = xt (1) 
where: 
t® — active usage time 
t™ — time at sea 
t — calendar time of the system observation 
«k=t™/t — time at sea factor (x € (0,1)) 


In view of these assumptions, the ICF type PS failures 
may occur only in the system active usage state, i.e. for the PS 
system in the t™ time, although their observed yearly numbers 
are determined by experts in relation to the calendar time t. 
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The model ICF probability has the vector form: 


i E Ak tE am, J 
Pit} =P{xt} = e, “k=l aK 0 
where: 

P{®} — the vector of probabilities of ICF type 
event occurrence within time interval 
(0,t) 


Ne > N; (t/t; — intensity function of HPP (ROCOF) 
a i (and at the same time the failure rate 

of the exponential distributions of time 

between failures in that process, [1/h] 


N, — annual number of the ICF type events 
elicited by j-th expert, [1/y] 

K; — time at sea factor elicited by j-th 
expert 

t; — calendar time of observation by j-th 
expert [h] 

J — number of experts 

K — the maximum number of possible ICF 


type failures in the time interval (0,t) 
t — the time of probability prediction. 


The A formula is based on the theorem on the asymptotic 
behaviour of the renewal process [1]: 


. ENG l 

iim eal 2] 2% (3) 
bone t To 

where: 

T, — mean time between failures. 


The number of ICF type events in the (0,t) period may be 
0,1,2,...or K with well-defined probabilities. The maximum of 
these probabilities is the assumed measure of the probability 
of ICF type event occurrence: 

(a),4\k a). 
Pmax{t} = max aro oe (4) 
ke{1,2,....K} k! 

The A and K parameters determined from the elicited 
opinions may be adjusted as new operation process data arrive 
on the investigated system failures. 

Expressions (2) and (4) allow to estimate the probabilities 
of ICF type hazardous events in the determined time interval t. 
Another problem is estimation of the risk of consequences 
of these events, i.e. damage to or total loss of the ship and 
connected human, environmental and financial losses. This is 
a separate problem not discussed in this paper. 


DATA ACQUISITION 


The PS will be further treated as a system consisting of 
subsystems and those consisting of the sets of devices. 

Experts are asked to treat the objects of their opinions 
as anthrop-technical systems, i.e. composed of technical 
and human (operators’ functions) elements. They elicit their 
opinions in three layers in such a way that proper correlation is 
maintained between data of the system and data of the system 
components. In layer 0 opinions are expressed in numbers, in 
layers I and II - in linguistic terms. For layers I and II separate 
linguistic variables (LV) and linguistic term-sets (LT-S) have 
been defined (Piegat A. 1999). 


Layer 0 — includes PS as a whole. 

Estimated are the annual numbers of type ICF type failures 
of PS N(t) and the percentage share of time at sea K100% in 
the time of expert’s observation. 


Layer I — includes decomposition of PS to a subsystem 

level. 

+ LV = share of the number of subsystem failures in the 
number of type ICF failures of PS. 

+ LT-S = Al-very small/none, Bl-small, Cl-medium, D1- 
large, El-very large. 


Layer II — includes the decomposition of subsystems to the 
sets of devices (set of devices is a part of subsystem forming 
a certain functional entity whose catastrophic failure causes 
catastrophic failure of the subsystem - e.g. it may be a set of 
pumps of the cooling fresh water subsystem). 

A LV = share of the number of failures of the sets of devices 
in the number of catastrophic failures of the respective PS 
subsystem. 

A LT-S = A2-very small/none, B2-small, C2-medium, D2- 
large, E2-very large. 


The structure of data acquisition procedure presented 
here implies a series form of the reliability structures of 
subsystems (layer I) and sets of devices (layer II). Elements 
of those structures should be so defined that their catastrophic 
failures cause equally catastrophic failures of the PS system 
and subsystem respectively. The division into subsystems and 
sets of devices should be complete and disjunctive. 

The data acquisition procedure presented here takes 
into account the expert potential abilities. It seems that their 
knowledge should be more precise in the case of a large 
operationally important system, as the PS is, and less precise 
as regards individual components of the system. 


ALGORITHM OF EXPERT OPINION 
PROCESSING 


In layer 0 the experts elicit annual numbers of the ICF type 
failures, which, in their opinion, might have occurred during 
1 year in the investigated PS type: 


N(t) j= 1,2,...53 (5) 


and shares of the time at sea in the calendar time of ship 
operation: 


«100% j = 1,2,...) 6) 
where: 
j — experts index 
J — number of experts. 


These sets of values are subjected to selection due to 
possible errors made by the experts. In this case a statistical 
test of the distance from the mean value may be useful, as in 
general we do not have at our disposal any objective field data 
to be treated as a reference set. 

If the data lot size after selection appears insufficient, it 
may be increased by the bootstrap method (Efron & Tibshirani 
1993). 

From the data (5) and (6), parameters A®and x of expression 
(2) and (4) are determined. Number of opinions J may be 
changed after the selection. 

In layer I experts elicit the linguistic values of subsystem 
shares in the number of ICF type failures of the investigated 
PS type (they choose LV value from the {A1, B1, C1, D1, El} 
set). The data are subjected to selection. 

The elicited data with linguistic values are compared in pairs 
- estimation of each subsystem is compared with estimation of 


each subsystem. The linguistic estimations are transformed into 
numerical estimations according to the following pattern: 
LT-S=Bl=>2 

LT-S=Cl>3 

LT-S=DI>4 

LT-S=El=>5 


Numerical estimates of each subsystem are subtracted 
from estimates of each subsystem. In this way the difference 
values are obtained, which may have the following values: 
-4,-3,-2,-1, 0, 1, 2, 3, 4. Those differences are transferred into 
preference estimates (as given in Tab. 1) in accordance with 
the following pattern: 

4 = 9, absolute preference 

3 = 7, clear preference 

2 > 5, significant preference 

1 = 3, weak preference 

0 => 1, equivalence 
-1 = 1/3, inverse of weak preference 
-2 > 1/5, inverse of significant preference 
-3 = 1/7, inverse of clear preference 
-4 > 1/9, inverse of absolute preference. 


From these differences, by the pair comparison method, 
a matrix of estimates is constructed. The estimates depend 
on the ,,distance” of the linguistic values LT-S of a given 
variable LV. For instance, preference A1 in relation to El has 
the value 9 assigned, in relation to D1 a value 7, in relation to 
C1 a value 5. In relation to B1 a value 3 and in relation to Al 
a value 1. The inverses of those preferences have the values, 
respectively: 1/9, 1/7, 1/5, 1/3 and 1. The matrix of estimates 
is approximated by the matrix of weight quotients of the 
sought arrangement. The recommended processing method is 
the logarithmic least squares method. The result is a vector of 
normalized arrangements of the subsystem shares (Saaty 1980, 
Kwiesielewicz 2002): 


P = [P;; Poss Piss Pil (7) 
where: 
p, — share of the i-th subsystem as a cause of an ICF type PS 
failure 


I — number of subsystems. 


Now we can determine in a simple way the intensity 
functions of individual subsystems arising from catastrophic 
failures: 


(a) _ 4(a) i= 
Ai =A" pi 1=1,2,...,1 (8) 
Tab. 1. Expert preference estimates acc. to Saaty (1980) 


Estimate Preference 


Equivalence 
Weak preference 
Significant preference 
Strong preference 
Absolute preference 


Inverse of the 
above numbers 


Inverse of the above described preference 


In layer II experts elicit the linguistic values of subset shares 
in the number of catastrophic subsystem failures (they choose 
LV value from the {A2, B2, C2, D2, E2} set). As in the case 
of subsystems, the expert opinions are processed to the form of 
normalized vectors of the arrangements of set shares: 


>The Saaty method, criticised in scientific circles, is widely applied in the decision-taking problems. 
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D-DaPet Peo Dl stole k= 129) 
where: 
p, — vector of the shares of i-th subsystem sets as causes of 


catastrophic failures of that subsystem 
share of the k-th set of i-th subsystem 
— number of sets in a given subsystem. 


Then, the intensity functions of sets contained in individual 
subsystems arising from catastrophic failures are determined: 


AW =A p. 3i=1,2,...,1;k =1,2,...,.K (10) 
EXAMPLE 


The example discusses investigation of a PS consisting 

of a low speed piston combustion engine driving a fix pitch 
propeller, installed in a container carrier ship. Experts were 
marine engineers with long experience (50 ship officers 
with chief engineer or second engineer diploma). Special 
questionnaire was prepared for them containing definition of 
the investigated object, schematic diagrams of subsystems and 
sets, precisely formulated questions and tables for answers. It 
was Clearly stated in the questionnaire that an ICF type failure 
may be caused by a device failures or by a crew actions. 
Out of 50 opinions elicited by experts, 3 were estimated as very 
unlikely (2 elicited numbers of the ICF events in a year were 
extremely underestimated and one was overestimated). They 
were eliminated and the remaining 47 opinions were further 
processed. 

Figs. 1 and 2 present statistical estimates of the expert 
opinion data (5) and (6). 
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Fig. 1. Box and whiskers plot of ICF yearly numbers 


Tab. 2 contains averaged basic data elicited by 47 experts 
in relation to the PS as a whole and the model parameters of 
ICF type event probability [equation (2)] determined from 
these data. 

From the Tab. 2 data the probabilities of determined numbers 
of ICF type event occurrences in | year were calculated. Fig. 3 
diagram presents results of those calculations. The numbers 
of probable ICF events in 1 year are equal 1, 2, ..., 5. The 
maximum probability is 0.2565, which stands for 2 ICF type 
events during | year, and the probability that such event will 
not occur amounts to 0.0821. 
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Fig. 2. Box and whiskers plot of time at sea share 


Tab. 2. Basic results of propulsion system investigation 


N(1 y) =2.5 
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Fig. 3. Disribution of ICF event numbers’ probability 


Tab. 3 contains the subsystem intensity function (ROCOF) 
data calculated from equation (8). The main PS risk 
“participants” are main engine and the electrical subsystem 
and the least meaningful is the propeller with shaft line. This 
is in agreement with the experience of each shipbuilder and 
marine engineer. 

Tab. 4 contains the fuel supply subsystem intensity function 
(ROCOF) data calculated from equation (10). 


Tab. 3. Intensity functions of the subsystems 


į 


Subsystem 


Fuel oil subsystem 


Sea water cooling subsystem 


Low temperature fresh water 
cooling subsystem 


Subsystem 2 10° 


High temperature fresh water 
cooling subsystem 


Starting air subsystem 


Lubrication oil subsystem 


Cylinder lubrication oil 
subsystem 


Electrical subsystem 


Main engine 


Remote control subsystem 


Propeller + shaft line 


(OPT 
a10 


Fuel oil service tanks 0.0488 | 2.2062 


Fuel oil supply pumps 7.5572 


Fuel oil circulating pumps 8.2840 


Fuel oil heaters. 4.2666 


Filters 6.9599 


10.6323 


Viscosity control arrangement 


Piping + heating up steam 


5.2965 
arrangement 


SUMMARY 


O The paper presents a method of subjective estimation of 


the hazard connected with losing by a seagoing ship of the 
propulsion function capability. The estimation is based on 
opinions elicited by experts - experienced marine engineers. 
The method is illustrated by an example of such estimation 
in the case of a propulsion system with a low speed diesel 
engine and a fix pitch propeller installed in a container 
carrier. 


O The given in section 6 do not raise any objections. The 


author does not have at his disposal sufficient objective data 
to evaluate precisely the adequacy of those data. It has to be 
taken into account that results of a subjective character may, 
by virtue of the fact, bear greater errors than the objective 
results achieved from investigations in real operational 
conditions. 

The presented method may be used in the procedures of 
the ship propulsion risk prediction. It allows to investigate 
the impact of the PS system component reliability on the 
probability values of ICF type event. It may also be used 
with other types of ship systems and not only to ship 
systems, particularly in the situations of hazardous event 
probability estimations with insufficient objective data at 
hand. 
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The application of the ahp method 
in ship system risk estimation 


Hoang Nguyen, Ph. D. 
Gdynia Maritime University 


ABSTRACT 


INTRODUCTION 


In sea transport, among numerous factors affecting the 
safety of navigation, of high importance is the reliability of the 
propulsion system. The situation when a ship losses its driving 
ability can lead, in some circumstances, to serious marine 
accidents, the consequences of which can be deaths of people 
and/or environment pollution. 

The risk of appearance of such an event, which we can 
refer to as the propulsion risk, depends on the reliability 
of individual components of the propulsion system and its 
operators. Estimating the risk of this system is troublesome due 
to the complexity of the system and the lack of historical data 
on its reliability. In such cases subjective estimations based on 
opinions of experts turn out very applicable [4]. 

In the article, the experts were ship engineers with many 
years’ experience who revealed numbers of propulsion system 
failures per year and linguistic proportions of contributions of 
particular subsystems to the total number of disastrous failures 
of this system. The estimation of risk function parameters for 
the subsystems included in the data presented by the experts 
was done by comparing in pairs. This way a proper correlation 
was obtained between the data referring to the entire system 
and its particular components (subsystems). The data obtained 
from comparing in pairs were processed using the Analytical 
Hierarchical Decision Making Process (AHP) method. 


THEORETICAL BACKGROUND 


The AHP method is a method commonly used in multi- 
attribute decision-making. The decision making process in 
the AHP method consists in ordering of variants, with further 
aggregation of the obtained ordered sets according to the 
defined hierarchical structure. The ordered sets are created via 
comparing in pairs. 
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The article presents procedures of application of the AHP method for estimating risks of 
ship systems. The estimation of system risk function parameters is done based on the data 
obtained from experts and processed using the pairwise comparative method. 


Key words: AHP method, estimation, propulsion risk, expert courts, ship systems 


The pairwise comparative method consists in comparing 
individual variants with each other and attributing a number 
from an earlier defined scale to each pair. The comparison is 
done by an expert who gives a precise (numerical) or fuzzy 
estimation to each pair. This estimation expresses individual 
expert’s preference when comparing one variant to the other 
(Tab. 1). 


Tab. 1. Estimations corresponding to experts’ preferences, 
after Saaty (1980) 


Estimation Preference 


Equivalence 


Weak preference of ith variant over 
jth variant 
Relatively strong preference of ith 
variant over jth variant 


Definite preference of ith variant 
over jth variant 
Absolute preference of ith variant 
over jth variant 


2,4, 6,8 Intermediate numbers 


Inverse numbers to 
those above listed 


Inverse preferences to those 
abovenamed 


Let us assume that we have n variants W,,W.,,...,;W, which 
are to be ordered. To each pair of variants (W, ,W,) an expert 
attributes a number r from the set S = {1/9, 1/8,. , 1/2, 1, 2, 

., 8, 9} which expresses his/her individual preferences with 
respect to the ith variant as compared to the jth variant. 


Experts’ estimations are used for creating the estimation 
matrices R (which are, as a rule, inconsistent). 


M D2 =- Nn 
r: > È 
Dhi In > 2n 
R= [I= : : : : (1) 
Th T2 =e Tan 


Matrix R is the matrix with inconsistent estimations and 
has the following aaa 


5,>0;%=+ 


> -ij 


Yi ;J= L2; 
ri 
Matrix R is considered consistent when its elements meet 
the condition: 


Kilik =Tik ; Vi, j,k = 1,2,...,n 


Ordering of variants is obtained as a result of approximation 
of the estimation matrix R using matrix P: 


Pu Pi +e Pin ) 


Pai P22 7 i (2) 


= = [p;]= T 
Pri Pn2 Parm 


the elements of which are consistent estimations presented in 
the form of weight ratios between particular variants: 


=i, 
Di: a Pad Ye 
where p, represent weights of particular variants of the order 
vector p: 
p= (Pp tee p3)” (3) 
After arithmetic standardisation we obtain the standardised 
order vector: 


p= (pi 2e pa (4) 


where: 


i=l 
Three basic methods which are used for approximating the 
estimation matrices are the following: 
* maximum eigenvalue method [3] 
* least square method [3] 
* logarithmic least square method [2, 3]. 


a) The maximum eigenvalue method 


This method consists in finding the vector p which fulfils 
the relation: 


Rp =\,,,,P (5) 


max: 
where: 
A — the maximum eigenvalue of matrix R. 


max 


b) The least square method 


This method consists in determining a form of matrix P 
which is closest to matrix R and which fulfils the relation: 


2 


min >. 5, —Pb (6) 


c) The logarithmic least square method 


This method consists in determining a form of matrix P 
which is closest to matrix R based on the Euclidean norm in 
the logarithmic scale. Matrix P fulfils the relation: 


2 
min; >?) Infr;)- In Pi (7) 
i,j= j 
After substituting: 
yi = In(r,) »x = In(p,) (8) 


we arrive at the optimisation problem: 
z 2 
minj I= $ [yj-x;+x;] 
i, j=l 


the solution of a has the form: 


4.25 a x +x;]= 0 (9) 
j k 
= i=l 
After introducing the condition of geometrical 
standardisation: 
n n 
[]pi=l; ox) =9 (10) 
i=l jl 
Equation (9) takes the form: 
n n 
1 f 
j=1 j=l 


Finally, after taking into account (8) we get the variant 
order vector: 
i n n 
j=l 


(12) 


For a large number of experts estimating a given pair 
of variants, the estimations made by individual experts are 
considered equivalent to each other. Then the task is limited 
to the following form: 


[ie 2 

min- +> In(y.)—In + (13) 
|i j=l k=1 j 

where: 

m= the number of estimations concerning the pair (i,j). 


After introducing substitutions: 


Vi, =In(r > x, = In(p,) (14) 
we arrive at the optimisation problem: 
n Mj 2 
i j=l k=1 
the solution of which has the form: 
d pea 
a“ 22 am +x,]=0;1=1,2,...n (15) 
at 


The set of equations (15) is transformed to the set of 
standard equations: 


n n n mij : 
x, ò m; -— X mix; = X X Yik E AA a 


j=l, j+i j=l, ji j=1,j+i k=1 


(16) 
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In matrix notation the set (16) has the form: 


Ax=b 
where: 
n mij 
b=[b;]; b;= >) Doviesi=L2...n 7 
j=1,j+i k=1 
>; —M)2 ~My 
j=l.j#l 
n 
A= 70i 2 mj =Ma (18) 
j=l. j#2 
n 
=s Ma -Ma2 ene $m 


j=l, jen 


Hence, the variant order vector can be determined from 
the equation: 


'b (19) 


[2] has shown that, despite its popularity, the maximum 
eigenvalue method reveals a weakness, namely it does not have 
the property of relation commutation, which means that the 
operation of estimation matrix transposition leads to different 
results. At the same time the least square method does not 
produce a unique solution. 

The logarithmic least square method is reduced to the 
geometric mean method for a large number of experts, as 
well as to the estimation matrix without missing data. The 
geometric mean method is a symmetrical transformation and 
produces the unique solution, irrelevant of scale inversion 
or changes in the order of particular aggregation operations. 
It is compatible with the estimation method of the highest 
credibility, which gives the reasons for its use from the point 
of view of the statistics. Further in the article, the logarithmic 
least square method is used for estimating the subjective ship 
propulsion risk. 

Among presently available computer codes which make 
use of the AHP method for decision-making, the codes which 
are preferred are Expert Choice and Criterium Decision Plus. 
Due to limitations concerning the numbers of variants in those 
codes, a computer code was developed based on calculating 
algorithms which make use of the logarithmic least square 
method. 


p=e* =e" 


APPLICATION 


The AHP method was used for estimating proportions 
of contributions of individual subsystems and units in the 
reliability structure of the ship propulsion system. The data 
were obtained from a group of 47 experts (ship engineers) using 
a specially prepared questionnaire [ESREL 2008]. 

The experts revealed their opinions on the proportions 
of contributions of particular subsystems in the total number 
of failures of the examined system in the form of linguistic 
values (very small, small, medium, big, very big). These values 
express experts’ preferences concerning the contributions 
of particular subsystems as the causes of propulsion system 
failures (Tab. 2). 

These data are compared in pairs and scaled using the 
assumed 5-stage scale to create the estimation matrices. 
A sample matrix which was created for the k-th expert is given 
below: 
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subsystems in total number of system failures. 


Name of 
subsystem 


Fuel oil 
subsystem 


very 
small 
/zero 


small 


medium 


Tab. 2. Linguistic experts’ opinions on the contributions of particular 


Seawater 
cooling 
subsystem 


Low 
temperature 
freshwater 
cooling 
subsystem 


High 
temperature 
freshwater 
cooling 
subsystem 


Starting air 
subsystem 


Lubrication 
oil 
subsystem 


Cylinder 
lubrication 
oil 
subsystem 


Electric 
power 
subsystem 


Main 
engine 


Main 
engine 
remote 
control 

subsystem 


.= [ty ]= 


Propeller 
and shaft 
line 


RK RANK KKK 


—_ 


æ Oo 


W = = — WT 
wr 


w 
Ww LW 


Nn nO = 
hn 


w 


-NX NX 


w uv N = 


— 


— 


— 


oS) 
an 


-X NXX 


Then the estimation matrices R are approximated using the order vector component ratio matrix. The estimation matrix 
approximation was performed using the logarithmic least square method for the case of a large number of experts. 
Matrix A is obtained in the form: 


l l l l l l 1 1 1 l 1) 
-47 470 -47 -47 -47 -47 -47 -47 -47 -47 -47 
-47 -47 470 =47 -41 -41 =47 -47 -47 =47 -47 
-47 -47 -47 470 -47 -471 -47 -47 -47 -41 -47 
-47 -47 -47 -47 470 -47 -47 -47 -47 -47 -47 

A=[a;]=|-47 -47 -47 -47 -47 470 -47 -47 -41 -41 -47 

47 -47 -47 -47 -47 -47 470 -47 -47 -47 -47 
-47 -47 -47 -47 -47 -47 -47 470 -47 -47 -47 

47 -47 -47 -47 -47 -47 -47 -47 470 -47 -47 
-47 -47 -47 -47 -47 -47 -471 -47 -47 470 -47 
-47 -47 -47 -471 -47 -47 -47 -47 -47 -47 470 


Vector b, vector x and the weight vector of standardised 
proportions of contributions of particular subsystems as causes 
of propulsion system failures, which were determined based on 
opinions obtained from a group of 47 experts, are as follows: 


0 0.5812 0.1329 

— 274.95 —0.5318 0.0436 
=327.13 — 0.6327 0.0394 
—94.04 —0.1819 0.0619 
71.12 0.1375 0.0853 
b=| -40.97 | ,x=| -0.0792 |, p* =| 0.0686 
— 264.76 -0.5121 0.0445 
478.40 0.9253 0.1876 
508.06 0.9827 0.1986 
212.73 0.4114 0.1122 
— 568.94 — 1.1004 0.0247 


Fig. 1 presents the estimated subsystem weights in the 
propulsion system function loss risk model, which were 
determined based on opinions given by 47 experts. 

It was assumed that the system unreliability function 
has exponential probabilities of the time to first failure. The 
exponential distribution of this time is characteristic for normal 
operation of numerous classes of systems, including ship 
systems (Gniedienko B.W et al., 1965). 


R(t) = Pfr <t)$=1-e; t>0 


where: 
À 
T = 


system failure rate 
random variable representing the time to failure. 


This model takes also into account standstills which are 
connected with renewals of the ship and its equipment and 
which are negligibly short compared to the ship maintenance 
time. Then the system failure rate À can be evaluated from the 
theorem about asymptotic behaviour of the recovery process 
(Gniedienko B.W. & Bielajew J.K. & Sotowiew A.D. 1965): 


No. of subsystem 


Fig. 1. Subsystem weights determined using the logarithmic least square 
method for the opinions obtained from 47 experts 


fig 
me t h 
where: 
T — mean time to failure 
N(t) — number of system failures. 


From the data revealed by the experts we get à = 3.34575 
E -— 04 t/h. 

The obtained subsystem contribution order vector can be 
written in the form: 


P=)" = (Pp Doce Pp aD) 


where: 

p, — contribution ofthe ith subsystem as the cause of system 
failure 

n — number of subsystems. 

The system failure rate can be determined from the 

relation: 


ai 


Lsi hi 
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The results of the calculations are shown in Tab. 3 as the CONCLUSIONS 
estimated weights and subsystem failure rates in the propulsion 
system function loss risk model. O The article presents the procedures of application of the 

pairwise comparative method to the ship system risk 


Tab. 3. Subsystem failure rates estimation. These methods are extremely applicable when 


Subsystem 


Fuel oil subsystem 


0.132970 


4.52026 


Seawater cooling 
subsystem 


0.043692 


1.48517 


Low temperature 


estimating expert data, which are, as a rule, inconsistent and/ 
or burdened with error. Although in general these methods 
are well known and in common use, the novelty of the here 
presented work consists in the application of the logarithmic 
least square method to ordering variants, which is reduced 
to the geometrical standardisation. The applied procedure 
secures proper correlation between opinions presented by 
many different experts. In case some data are missing this 
method can also be used directly, and the missing data are 
assumed as consistent with the resultant order. 


freshwater cooling 0.039497 1.34260 O Currently, activities are in progress in the Department of 
subsystem Engineering Sciences, Gdynia Maritime University, upon 
ore ae ea a neuron-fuzzy model of the ship propulsion risk. The 
fre a ae 0.061996 2.10739 estimated parameters obtained from the expert base using 
subsystem 8 i , the AHP method are the input data for this model. 
BIBLIOGRAPHY 


Starting air 
subsystem 


0.085332 


2.90062 


Lubrication oil 
subsystem 


0.068698 


2.33519 


Cylinder lubrication 
oil subsystem 


0.044561 


1.51474 


Electric power 
subsystem 


0.187603 


6.37703 


Main engine 


Main engine remote 
control subsystem 


0.198681 


0.112219 


6.75360 


3.81457 
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Propeller and shaft 
line 
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Supporting navigator’s decisions 
by visualizing ship collision risk 
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ABSTRACT 


The paper introduces a visualization method that enables the navigator to estimate an encounter situation 
and choose collision avoidance manoeuvre, if necessary. It is based on the Collision Threat Parameters 
Area method and offers new features: fuzzy sectors of forbidden speed and course values, the possibility 
to use any given ship domain and a new formula of collision risk assessment. All these elements result in 
a method, which enables the navigator to differentiate between varying levels of risk and to point out the 
direct threats. The method is fast enough to be applied in the real-time decision support system. 


Key words: ARPA display, collision avoidance manoeuvres, collision risk, navigation 


INTRODUCTION 


Traditional displays used in collision-avoidance systems 
were based on the relative Cartesian coordinate system, with 
the own ship in the centre of it and X and Y coordinates 
denoting the relative positions of all other objects including 
ships — here called targets. Their functionality was limited 
to showing all targets within a certain range and indicating 
the targets that were considered to be dangerous on the basis 
of computations performed by the system. Some of them 
additionally visualized Predicted Areas of Danger (PAD) 
and the resulting necessary course alterations [9]. What these 
displays failed to visualize was the nature of collision risk: 
the colliding combinations of courses and speeds of the own 
ship and the dangerous targets. Visualizing these forbidden 
combinations of course and speed (instead of course only) 
has been introduced by Lenart as a part of Collision Threat 
Parameters Area (CTPA) method [1]. 

Based on CTPA the graphical method determining necessary 
manoeuvre has been introduced in [6, 7, 8]. However, the display 
according to Lenart and the derived method naturally assume 
using a pre-defined safe distance D, as a main safety parameter 
and consequently DCPA (Distance to the Closest Point of 
Approach) as a collision risk measure. Therefore they cannot 
be used for precise visualization of the necessary manoeuvres 
when domains other then circle-shaped are assumed. Also, in 
cases of large values of safe distance parameter D, and multi 
ship encounters it might be impossible to find one combination 
of course and speed that avoids collisions with all targets. In 
such cases the targets should be ranked according to their risk 
factor values and a sequence of collision avoidance manoeuvres 
should be applied. 


In [5] the following ideas have been combined: double 
coordinate system used in CTPA, a fuzzy ship domain [3] and 
approach factor f „~ a new measure of collision risk derived 
from the concept of a ship domain [4]. Approach factor f „has 
been defined as the scale factor of the largest domain-shaped 
area that is predicted to remain free of other ships throughout 
the whole encounter situation. The result presented in [5] has 
been a new visualization method called Fuzzy Collision Threat 
Parameters Area (FCTPA). It has been based on Collision 
Threat Parameters Areas method and has extended it so as to 
handle any given domain, including fuzzy domains. 

In this paper the FCTPA method is further modified and 
a generalized version of a formula for collision risk [6, 7, 8] 
is applied. The original formula presented in [6, 7, 8] took 
into account both DCPA (Distance at the Closest Point of 
Approach), and TCPA (Time remaining to the Closest Point 
of Approach), thus emphasizing the direct threat of the nearest 
close encounters. It also reflected the non-linearity of the 
collision risk (as a function of distance and time) near the 
critical values of minimum safe distance D, and minimum time 
for manoeuvring T,. The generalized version of the formula 
replaces DCPA and TCPA with f „ and t — the approach 
factor value and the time remaining to reaching this value. 
The modified version of the FCTPA visualization method 
applying the new formula for collision risk factor is called 
Two-Parameter Collision Threat Area (TPCTA). 

The rest of the paper is organized as follows. Section 2 
briefly describes the CTPA method. In Section 3 the formulas 
for collision risk assessment are provided and in Section 4 the 
new TPCTA method is presented. An exemplary scenario of 
using this method is given in Section 5. Finally Section 6 is 
a summary of the paper. 
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COLLISION THREAT PARAMETERS 
AREA (CTPA) 


In [1] a collision threat is defined as a target ship for which 
the following condition holds: 


DCPA < D, (1) 


The method uses a double Cartesian coordinate system 
where the horizontal axis represents both the X coordinate 
of position and V, coordinate of speed and the vertical axis 
represents both the Y coordinate of position and V, coordinate 
of speed. The relation between the position and speed 
coordinates is as follows: 


x=V *T 
(2) 
<= Vo a 
y 
where: 
t — fixed time value, for example 12 minutes. 


The Collision Threat Parameters Area for a single target 
ship is defined as an area in the abovementioned system of 
coordinates that fulfils the following conditions: 

* placing the tip of the own speed vector V within this area 
would result in violating the safe distance D, between 
the ships 

* placing the tip of the own speed vector V outside this 
area would result in keeping the safe distance D, between 
the ships. 


The Collision Threat Parameters Area for a group of target 
ships is defined as a superposition of the Collision Threat 
Parameters Areas obtained for each of the targets separately. 
Graphical superposition has been presented in [6, 7, 8]. 
The formula for the two straight lines determining the boundaries 
of the CTPA for a given single target is as follows: 


y=ax—b,1 


(3) 
y=ax=—b,t 
where the coefficients are given by the formulas: 
(4) 
(5) 
where: 
x,,y, — coordinates of the relative position of the target 


ship 
Vas Vos coordinates of the true speed of the target ship. 


In practice, CTPA is only this part of the determined 
area, where the condition TCPA > 0 holds, since only future 
collision threats are of interest. Also, in case of a multiple 
target encounter, the manoeuvres for which the safe distance D, 
would be violated after a time longer than the critical time 
(DCPA < D,, TCPA > T..) may be allowed, if there is no 
possibility of avoiding all targets with just one manoeuvre. 
This means that the tip of the own speed vector may be 


84 POLISH MARITIME RESEARCH, No 1/2009 


conditionally placed within this part of the CTPA, for which 
TCPA > T.. When applied to the graphical display, the CTPA 
method enables the operator to choose manually a safe own 
speed vector in a very easy way — it is enough to choose a point 
outside the CTPA and read its speed coordinates. The method 
is summarized by Fig. | and Fig. 2. 


K 


Fig. 2. Graphical interpretation of an encounter of two ships: 
before a collision avoidance manoeuvre (up side) 
and after the manoeuvre (down side) 


The line of DCPA = 0 represents a perfect crash. If the tip 
of the own speed vector is placed above this line, it means 
that the own ship will pass ahead of a target, if below this line 


— astern of a target. Fig. 2 presents a graphical interpretation 
of an encounter of two ships: before and after a collision 
avoidance manoeuvre. The own ship is moving with the speed 
of 15 knots and course 0 degrees and a target ship with the 
speed of 10 knots and course 300 degrees. The tip of the own 
speed vector is placed above the line for DCPA= 1.5 n. m. 
Therefore a collision avoidance manoeuvre is taken and the 
new own course is set to 35 degrees. New positions of ships are 
visualised after the time t representing the physical realization 
of this manoeuvre. 


POINTING OUT THE DIRECT THREATS IN 
MULTI SHIP ENCOUNTER SITUATIONS 


To determine which targets constitute direct threat and 
have to be avoided first, it is necessary to monitor their motion 
parameters and to assign to them values of collision risk factor 
according to a given formula. Such formula for collision risk 
factor has been proposed by Lisowski [2] and it is as follows. 

l 


pcpA)  (TCPA ) DYI? 
r=|a, +a, +a] — (6) 
Ds i Ty l Ds 


collision risk factor 


D — current distance between the own ship and the 
target ship 

Di — safe distance of approach (a radius of the circle- 
shaped domain), 

Ts — asafe time necessary to plan and perform a collision 
avoidance manoeuvre, 

a,,a,,a, — weight coefficients, dependent on the state of 


visibility at sea, dynamic length and dynamic beam 
of the ship and a kind of water region. 


The advantage of this formula is that it can be applied to 
any target. However, the formula will indicate a considerable 
risk for DCPA = D, and TCPA = 0, that is, when the ships pass 
each other at a safe distance and there is no risk involved any 
more. 

It must be emphasized that the most important elements 
affecting collision risk are: vector of relative speed (the 
superposition of the own speed and the speed of a target) and 
the lack of change of the target’s bearing. These elements have 
been taken into account in [6,7,8]. Additionally it has been 
proposed there to limit the considered range of DCPA and 
TCPA values to the real collision risk area.(| DCPA | € (0, D) 
and TCPA € (0, nT .), as given by conditions of formula (7). 
CTPA method by Lenart has been used, with the following 
modification applied: the function (7) has been introduced, 
which is strongly non-linear around the critical values of DCPA 
and TCPA. Thus modified collision risk factor r as a function 
of DCPA and TCPA is also shown in Fig. 3. 


=, T, DCPA < D 
‘ y “PA < D, 
a*le | Ps l So J for ad 
TCPA TCPA < nT, 
(7) 
r=; DCPA > D; 
| 0 for or i 
TCPA 2 nT, 
where: 
1>r>0n>1 
a,b,c — values determined in the course of simulation: 


a=1.11,b=1.52,c=0.33. 


Fig. 3. Geometrical interpretation of collision risk factor 


In the formula above the DCPA/D, quotient is used, which 
is exactly the approach factor value f „ [4] for a circle-shaped 
domain. In case of using a ship domain model other than that, 
this quotient can be replaced with the f „„ value obtained for this 
model. Similarly TCPA might be replaced with the time remaining 
to reaching f ,, value. The generalization of the formula (7), which 
would hold true for any given ship domain would thus be: 


Ts fain < l 


pe. bE 
a*le fmin 0,1] —c| for 


fmin 
r=; f min 2 l 
lo for or 
TE ain 2 nTg 
Using the formula (8), we can rank the targets according to 
the decreasing values of the risk factor r. 


DIRECT COLLISION THREAT 
PARAMETERS AREA (DCTPA). 


In this section a visualization tool that has been designed by 
the authors — Direct Collision Threat Parameters Area (DCTPA) 
is presented. It is based on the same concept of the forbidden 
area in the double Cartesian coordinate system, but instead 
of determining the area analytically for a fixed D, value, it is 
determined numerically for a given ship domain. It works as 
follows. For every combination of the own course and speed 
the resulting f „ value is computed (the algorithms used to 
compute the f value for given courses, positions and speeds 
of the own ship and target ships have been presented in detail 
in [4]). Then for each combination of the own course and speed 
the risk factor for each target is computed according to formula 
(8) and the minimal of those risk factor values is assigned to 
the current combination of course and speed. 

Depending on the obtained value, a point in the double 
Cartesian coordinate system representing a particular 
combination of the own speed and course is assigned a colour 
in the following way: 
> for 0.75 <r <1 (critical collision risk): black 
> for 0 <r < 0.75 (significant collision risk): gradually 

changing dark grey 
> for 0 < r < 0.5 (significant collision risk): gradually 

changing light grey 
> for r=0 (safe passing): white 
The proposed visualisation tool is shown in Fig. 4 
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Fig. 4. Direct Collision Threat Parameters Area 


Whenever the arrow indicating the end of the own 
speed vector appears on the dark grey or black background 
— a collision avoidance manoeuvre should be performed. The 
closest white point on the display represents a safe combination 
of the own speed and course. For the situation in Fig. 4 the tip 
of the own speed vector is on the white background (outside the 
forbidden sector associated with the target), hence the current 
combination of course and speed is safe. 


AN EXEMPLARY SCENARIO 


In this section it is shown, how visualizing collision risk 
affects navigational decision. An example, randomly generated 
scenario, involving 20 targets is analysed here. The initial 
situation is given in Tab. | and Fig. 5. 

According to Tab. 1 and Fig. 5, the own ship is approaching 
a number of target ships, some of them colliding with the 


Tab. 1. The positions, courses and speeds of the own ship and the target ships in the exemplary scenario 


Course 
[degrees] 


Position coordinates at the start 


Coordinates of the 


time destination point 


Own ship 
Target 1 


x [n.m.] 


y [n.m.] 


Target 2 


Target 3 


Target 4 


Target 5 


Target 6 
Target 7 


Target 8 


Target 9 


Target 10 


Target 11 


Target 12 
Target 13 


Target 14 


Target 15 


Target 16 


Target 17 


Target 18 
Target 19 


Target 20 
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Fig. 5. Position and courses of the own ship and target ships in the exemplary scenario 
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Fig. 7. Forbidden sectors of own course and speed for the exemplary scenario according to DCTPA 
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Fig. 8. A safe trajectory of the own ship resulting from a manoeuvre to starboard 


own course. The forbidden sectors of own course and speed 
(Collision Threat Parameters Area) for the safe distance of 0.5 
nautical mile look as shown in Fig. 6. Due to the large number 
of targets it is impossible to choose a single course alteration 
manoeuvre to starboard (recommended by COLREGS) that 
avoids collisions with all of the targets. 

However, applying the proposed visualization tool enables 
the navigator to assess the situation better than that. According 
to the proposed method, the Direct Collision Threat Parameters 
Area for the initial situation looks as depicted In Fig. 7. 
The Coldwell domain is used here for target ships. 

From the collision risk visualized in Fig. 7, the navigator 
may easily conclude, that although the current course and speed 
are not safe (the tip of the own speed vector is on the forbidden 
area), itis possible to perform a course alteration manoeuvre to 
starboard, which avoids collision with the directly threatening 
target. As a result, a trajectory containing the manoeuvre to 
starboard might be planned, as shown in Fig. 8. 


SUMMARY AND CONCLUSIONS 


Visualizing collision risk is especially important in 
multi-ship encounter situations, when the own course and 
speed collides with more than one target. Unfortunately, in 
such situations it is not always possible to avoid collision by 
performing a single, lawful collision avoidance manoeuvre. 
Therefore the targets have to be ranked according to the risk 
they constitute. This risk depends not only on the minimum 
distance between the ships (or the domain violation) but also 
on the time remaining to the point of the closest approach. 
Particularly threatening are situations when domain violations 
are accompanied by short times remaining to the encounters. In 
such situations the proposed solution enables the navigator to 
point out the most dangerous targets and avoid them in the first 
place. Once these targets have been safely passed, the navigator 
may focus on the targets constituting lesser risk. The proposed 
solution is fast enough to be applied in the real-time decision- 
support system, where fast processing of the data concerning 
all targets is a necessity. 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


Ə Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 


Iława Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 


The Foundation possesses ISO 9001 quality certificate. 
Why training on ship handling? 


The safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training is to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office@ilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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GDANSK UNIVERSITY OF TECHNOLOGY 


is the oldest and largest scientific and technological academic institution in 
the Pomeranian region. The history of Gdansk University of Technology is 
marked by two basic dates, namely: October 6, 1904 and May 24, 1945. 


The first date is connected with the beginning of the technical education at 
academic level in Gdansk. The second date is connected with establishing of 
Gdansk University of Technology, Polish state academic university. Gdansk 
University of Technology employ 2,500 staff, 1,200 whom are academics. 
The number of students approximates 20,000, most of them studying 


full-time. Their career choices vary from Architecture to Business and (© O a f) S Ki U f) iva IS jy 
of Tacnnology 


Management, from Mathematics and Computer Science to Biotechnology 
and Environmental Engineering, from Applied Chemistry to Geodesics and 
Transport, from Ocean Engineering to Mechanical Engineering and Ship 
Technology, from Civil Engineering to Telecommunication, Electrical and 
Control Engineering. Their life goals, however, are much the same - to 
meet the challenge of the changing world. The educational opportunities dia 
offered by our faculties are much wider than those of other Polish Technical + 
universities, and the scientific research areas include all of 21st Century F It f @) woio 

technology. We are one of the best schools in Poland and one of the best d G U © G € d f >} 
known schools in Europe — one that educates specialists excelling in the 


programming technology and computer methods used in solving complicated 4 4 PG 
scientific, engineering, organizational and economic problems. E f g | f 5 H id f g we 
THE FACULTY OF OCEAN ENGINEERING 
AND SHIP TECHNOLOGY 


The Faculty of Ocean Engineering and Ship Technology (FOEST) as 
the only faculty in Poland since the beginning of 1945 has continuously 
been educating engineers and doctors in the field of Naval Architecture 

and Marine Technology. 


The educational and training activities of FOEST are supported by 
cooperation with Polish and foreign universities, membership in different 
international organizations and associations, as well as participation in 
scientific conferences and symposia. Hosting young scientists and students 

from different countries is also a usual practice in FOEST. 


The activities of Faculty departments are related to: mechanics and 
strength of structures, hydromechanics, manufacturing, materials and system 
quality, power plants, equipment and systems of automatic control, mostly 

in shipbuilding, marine engineering and energetic systems. 


FOEST is a member of such organizations like WEGEMT; The 
Association of Polish Maritime Industries and the co-operation between 
Nordic Maritime Universities and Det Norske Veritas. The intensive teaching 
is complemented and supported by extensive research activities, the core 
of which is performed in close collaboration between FOEST staff and 
industry. We take great care to ensure that the applied research meet both 
the long term and short term needs of Polish maritime industry. FOEST 
collaborates with almost all Polish shipyards. Close links are maintained 
with other research organizations and research institutions supporting the 
Polish maritime industry, such as Ship Design and Research Centre and 
Polish Register of Shipping, where several members of the Faculty are also 

members of the Technical Board. 


The Faculty of Ocean Engineering and Ship Technology is a unique 
academic structure, which possesses numerous highly qualified and 
experienced staff in all above mentioned specific research areas. Moreover, 
the staff is used to effective co-operation and exchange of ideas between 
specialists of different detailed areas. This enables a more integrated and 
comprehensive treatment of research and practical problems encountered 
in such a complicated field of activity as naval architecture, shipbuilding 

and marine engineering. 


The staff of the Faculty has strong international links worldwide, being 
members or cooperating with international organizations like International 
Maritime Organization IMO, International Towing Tank Conference ITTC, 
International Ship and Offshore Structures Congress ISSC, International 
Conference on Practical Design of Ship and other floating Structures PRADS 

just to name a few. 
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